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Abstract 

The  purpose  of  this  research  effort  is  to  develop,  simulate,  and  test  a  new 
algorithm  to  detect  Near  Earth  Objects  (NEOs)  using  a  Likelihood  Ratio  Test  (LRT) 
based  on  a  Poisson  statistical  model  for  the  arrival  of  photons.  One  detection  algorithm 
currently  in  use  is  based  on  a  Gaussian  approximation  of  the  arrival  of  photons,  and  is 
compared  to  the  proposed  Poisson  model.  The  research  includes  three  key  components. 
The  first  is  a  quantitative  analysis  of  the  performance  of  both  algorithms.  The  second  is  a 
system  model  for  simulating  detection  statistics.  The  last  component  is  a  collection  of 
measured  data  to  apply  comparatively  to  both  algorithms. 

A  Congressional  mandate  directs  NASA  and  the  DoD  to  catalogue  90%  of  all 
NEOs  by  the  year  2020  [1].  Results  from  this  research  effort  could  feasibly  be  applied 
directly  to  operations  in  the  Pan-Starrs  program  to  facilitate  the  accomplishment  of  the 
Congressional  mandate.  Improvements  in  the  size  of  detectable  NEOs  and  in  the 
probability  of  detecting  larger  NEOs  would  increase  the  state  of  readiness  of  the  world 
for  possible  catastrophic  impact  events.  Improvements  in  detection  probability  of 
measured  data  were  as  high  as  a  factor  of  seven,  and  the  expected  average  improvement 
is  10%. 
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NEAR  EARTH  OBJECT  DETECTION  USING  POISSON  STATISTICAL 


MODEL  FOR  DETECTION  ON  IMAGES  FROM  THE  PANORAMIC  SURVEY 
TELESCOPE  &  RAPID  RESPONSE  SYSTEM 


1  Introduction 


1.1  General  Issue 

In  the  National  Aeronautics  and  Space  Administration  (NASA)  Multiyear 
Authorization  Act  of  1990,  the  United  States  Congress  directed  a  workshop  study  to 
define  a  program  to  increase  the  detection  rate  of  asteroids  whose  trajectory  crosses  the 
orbital  path  of  Earth.  This  lead  to  NAS  As  Spaceguard  Survey  Report  in  1992,  the 
conclusions  of  which  called  for  a  worldwide  network  of  4  to  7  telescopes  in  the  2  to  3 
meter  aperture  range.  The  report  predicted  that  nearly  all  asteroids  and  comets  over  1  km 
in  diameter  could  be  catalogued  and  tracked  with  such  a  network.  The  report  predicts 
that  10  percent  of  smaller  asteroids  and  comets  between  100  meters  and  1  km  in  diameter 
could  be  catalogued  and  tracked  with  a  similar  system.  The  report  outlines  the  need  for 
cooperation  with  the  Department  of  Defense,  particularly  with  the  US  Air  Force,  in  the 
search  for  NEOs  [2], 

In  1994,  the  House  Committee  on  Science  and  Technology  directed  NASA,  in 
coordination  with  the  Department  of  Defense  and  other  international  space  agencies  to 
discover,  catalogue,  and  track  within  10  years,  90  percent  of  all  asteroids  and  comets 
larger  than  1  km  within  1.3  Astronomical  Units  (AU)  from  the  sun  whose  trajectory 
crosses  the  orbital  path  of  Earth.  Eleven  years  later,  the  NASA  Authorization  Act  of 
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2005  issued  a  new  mandate  in  which  90%  of  all  NEOs  larger  than  140  meters  in  diameter 
must  be  discovered,  catalogued,  and  tracked  by  the  year  2020.  Congress  also  directed 
NASA  to  submit  an  Analysis  of  Alternatives  (AoA)  to  the  Committee  within  120  days  of 
enactment  of  the  Act,  outlining  efforts  taken  by  NASA  to  detect  and  characterize  the 
hazards  of  NEOs,  as  well  as  an  assessment  of  necessary  actions  to  put  in  place 
capabilities  to  expand  detection  and  tracking  of  NEOs  [1]. 

The  AoA  submitted  to  Congress  by  NASA  in  2007  details  two  considered 
terrestrial  detectors,  and  several  space-based  systems.  The  two  terrestrial  based  systems 
are  the  Large  Synoptic  Survey  Telescope  (LSST),  and  Panoramic  Survey  Telescope  and 
Rapid  Response  System  (Pan-Starrs)  [3],  The  AoA  reported  that  a  program  consisting  of 
a  combination  of  both  ground  based  systems  and  some  space-based  systems  was  required 
to  meet  the  2020  deadline  for  completion.  The  AoA  also  reported  that  using  only  one  of 
the  land-based  systems  would  push  the  date  out  to  beyond  2030.  The  number  of  NEOs 
estimated  by  the  AoA  is  depicted  in  Figure  1,  using  a  constant  power  law  [3]. 
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Figure  1.  Frequency  of  NEOs  by  Size,  Impact  Energy,  and  Magnitude  [3] 
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Figure  2.  Impact  Destruction  Radius  [4] 

Large  NEOs  have  an  impact  frequency  of  one  every  half  a  million  years,  while 
sub-kilometer  NEOs  impact  one  in  every  thousand  years  [4],  While  the  probability  of  an 
impact  of  a  large  NEO  capable  of  triggering  mass  extinctions  within  a  lifetime  is  small, 
the  probability  of  an  impact  from  a  smaller  NEO  is  significantly  higher,  and  as  Figure  2 
illustrates,  such  an  impact  would  cause  catastrophic  localized  damage. 

1.2  Problem  Statement 

NASA's  budget  request  for  NEO  observations  rose  from  $5.8M  in  FY  2010  to 
$20.4M  in  FY  2011  and  beyond,  the  largest  increase  since  NASA  began  searching  for 
NEOs  [5].  This  increase  reflects  the  realization  that  enough  is  not  being  done  to 
accomplish  the  Congressional  Mandate  outlined  in  [1],  The  cost  for  the  spaced-based 
systems  outlined  in  [3]  was  forecasted  to  be  more  than  twice  that  of  the  terrestrial -based 
systems,  as  well  as  being  more  complicated  to  support  and  maintain.  The  cost  of  the 
terrestrial-based  systems  was  forecasted  to  be  $469M  and  would  not  meet  the 
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Figure  3.  Total  NEO  Discoveries  [6] 

requirements  of  the  mandate  until  the  year  2026  [3],  As  of  7  October  2011,  fewer  than 
900  NEOs  larger  than  1  km  have  been  discovered,  and  fewer  than  5000  between  100  m 
and  1  km  have  been  discovered,  as  shown  in  Figure  3  [6],  With  the  given  operating 
budget  for  NEO  detection,  the  required  budget,  and  the  current  progress  of  the  effort  to 
complete  the  Congressional  mandate,  a  solution  is  needed  that  improves  the  detection 
capability  of  current  equipment  at  minimal  costs. 


1.3  Research  Objectives/Focus 

The  objectives  of  this  research  effort  will  be  to  develop  a  new  algorithm  for 
detecting  NEOs  using  existing  hardware,  namely  Pan-Starrs.  The  research  will  be 
focused  on  decreasing  the  detectable  size  and  increasing  the  probability  of  detecting 
larger  NEOs  by  changing  the  post-processing  algorithm  of  image  data.  The  model  for  the 
simulations  will  be  greatly  simplified  from  the  complex  algorithm  that  is  currently  in  use. 
All  of  the  functionality  of  Pan-Starrs  is  not  investigated  in  this  report;  therefore,  only  the 
process  by  which  Pan-Starrs  takes  image  data  and  uses  it  to  flag  the  detection  of  a  NEO  is 
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modeled.  Likewise,  the  investigated  algorithm  would  likely  require  additional 
functionality  and  integration  work  in  order  to  implement. 

1.4  Investigative  Questions 

How  does  Pan-Starrs  currently  detect  NEOs?  What  are  the  optical  capabilities  of 
Pan-Starrs?  Given  the  capabilities  of  the  Pan-Starrs  hardware,  what  detection  theory  can 
be  applied  to  the  post-processing  of  the  image  data?  How  does  any  proposed  detection 
algorithm  compare  to  the  existing  detection  algorithm?  Can  any  proposed  algorithm 
extend  to  applications  other  than  Pan-Starrs? 

1.5  Methodology 

Receiver  Operating  Characteristic  (ROC)  curves  will  be  generated  through  Monte 
Carlo  simulations  based  on  models  of  the  Pan-Starrs  optical  characteristics,  atmospheric 
characteristics,  and  NEO  characteristics.  Simulations  of  varying  environments  such  as 
proximity  to  brighter  objects,  and  the  intensity  of  nearby  objects  will  be  investigated. 
Tests  in  relative  environments  will  be  conducted  to  verify  the  results  of  the  simulations. 

1.6  Assumptions/Limitations 

Many  variables  are  involved  in  the  ability  to  detect  a  NEO.  Without  actual  data 
from  Pan-Starrs,  assumptions  must  be  made  about  the  environment  in  which  images  will 
be  taken.  Average  atmospheric  conditions,  average  NEO  characteristics,  and  the  average 
optical  response  of  Pan-Starrs  will  be  used  in  the  simulations.  Similarly,  without  the 
availability  of  Pan-Starrs  hardware,  a  relative  test  environment  will  be  derived  based  on 
the  results  seen  in  simulations. 
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1.7  Implications 

Results  from  this  research  effort  could  feasibly  be  applied  directly  to  operations  at 
Pan-Starrs  and  other  programs.  If  improvements  in  the  size  of  detectable  NEOs  or  in  the 
probability  of  detecting  larger  NEOs  are  observed,  the  algorithm  could  be  used  to 
facilitate  the  accomplishment  of  the  Congressional  mandate  to  have  90%  of  all  NEOs 
over  140  m  in  diameter  catalogued  by  2020.  Any  such  notable  improvement  would 
increase  the  state  of  readiness  of  the  world  for  possible  catastrophic  impact  events. 

1.8  Preview 

This  research  aims  to  demonstrate  that  detection  theory  can  be  used  to  implement 
a  Likelihood  Ratio  Test  (LRT)  to  improve  upon  the  current  detection  algorithm  used  by 
Pan-Starrs.  Further,  it  aims  to  demonstrate  that  other  known  electronic  filtering 
techniques  are  based  on  approximations  about  the  stochastic  nature  of  photon-counting 
and  that  improved  results  are  possible  without  making  such  approximations.  Chapter  2, 
Literature  Review,  compares  the  previous,  current  and  planned  NEO  detection  methods 
of  multiple  NEO  detection  programs.  Chapter  3,  Methodology,  details  the  analytic 
process  that  led  to  the  proposed  algorithm,  as  well  as  the  approach  used  to  develop 
models,  simulations,  and  tests.  Chapter  4,  Analysis  and  Results,  interprets  the  outcomes 
of  the  analysis,  simulations,  and  tests.  Chapter  5,  Conclusions  and  Recommendations, 
discusses  the  validity  and  perfonnance  of  the  proposed  algorithm  based  on  comparisons 
of  the  analytical,  simulation  and  test  outcomes  with  the  performance  of  current  NEO 
detection  methods. 
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2  Literature  Review 


2.1  Chapter  Overview 

The  purpose  of  this  chapter  is  to  review  past,  current,  and  future  NEO  detection 
methods.  Beginning  in  1984  with  the  Spacewatch  program,  through  the  estimated  2020 
operational  date  of  the  Large  Synoptic  Survey  Telescope  (LSST)  program  there  have 
been  numerous  government  and  academic  programs  dedicated,  at  least  in  part,  to  the 
discovery  and  cataloguing  of  NEOs.  An  analysis  of  the  algorithms  and  system 
capabilities  of  such  programs  will  show  that  the  proposed  algorithm  in  this  research  effort 
is  a  novel  approach  to  the  post-processing  of  astronomical  image  data.  NEO  detection 
programs  have  employed  increasingly  more  capable  and  sophisticated  optical  systems 
and  computer  processing  systems;  however,  very  little  has  changed  in  the  basic  method 
of  detecting  NEOs. 

2.2  Previous  and  Current  Programs 

In  1984,  University  of  Arizona’s  Spacewatch  program  became  the  first  to  use 
Charge-coupled  Devices  (CCDs)  to  scan  the  sky  for  NEOs.  The  Spacewatch  team 
developed  the  tools  and  methods  used  by  many  of  the  astronomical  observing  telescopes 
today  to  automatically  detect  moving  objects  using  CCDs  [7],  The  “drift-scan”  method 
focuses  the  telescope  at  a  point  slightly  leading  the  targeted  area  and  keeps  it  stationary, 
allowing  the  sky  to  drift  across  the  field  of  view  due  to  the  rotation  of  the  Earth.  The  read 
out  rate  of  the  CCD  is  set  equal  to  the  drift  rate  of  the  sky,  providing  a  very  long  exposure 
time  while  keeping  any  fixed  objects  in  focus.  Objects  that  are  moving  fast  relative  to  the 
background  sky,  as  NEOs  would  be,  would  appear  as  streaks  in  the  image.  The  “step- 
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stare”  method  uses  a  mosaic  of  fast  read  out  CCDs  to  piece  together  an  image  then  return 
to  the  same  area  of  sky  a  number  of  times  within  a  small  timeframe  in  order  to  perform 
some  type  of  change  detection  between  the  two  images.  NEOs  moving  faster  than  the 
background  sky,  instead  of  appearing  as  streaks,  will  appear  in  a  different  location 
relative  to  the  surrounding  objects  in  the  image  [8]. 

Since  the  Spacewatch  program,  at  least  five  subsequent  NEO  detection  programs 
have  used,  or  are  using  the  "step-stare"  method  in  order  to  detect  moving  objects.  These 
programs  include  the  Near  Earth  Asteroid  Tracking  (NEAT),  Lincoln  Near  Earth 
Asteroid  Research  (Linear),  Lowell  Observatory  NEO  Search  (LONEOS),  Catalina  Sky 
Survey  (CSS),  and  the  Japanese  Spaceguard  Association  (JSGA)  programs.  Each  of 
these  programs  has  offered  unique  capabilities  and  contributions,  but  a  detailed  review  of 
each  program  is  outside  of  the  scope  of  this  document  [9]  Table  1  displays  a  summary 
comparison  of  the  six  program  telescopes  previously  discussed  [8]. 


Table  1.  Comparison  of  current  NEO  search  programs 


Program 

Space- 

watchl 

Space- 

watch2 

NEAT-1 

NEAT-2 

LONEOS 

LINEAR 

Catalina 

Sky 

Survey 

JSGA 

Aperture 

(m) 

0.9 

1.8 

1.2 

1.2 

0.6 

1 

0.4 

0.5 

f# 

5.3 

2.7 

1.9 

2.5 

1.9 

2.2 

3 

1.9 

Pixel  Size 
(mm) 

0.024 

0.024 

0.015 

0.015 

0.0135 

0.024 

0.015 

0.015 

Pixel  Size 
(arcsec) 

1 

1 

1.4 

1 

2.5 

2.25 

2.5 

3.2 

FOV  (deg2) 

0.3 

0.3 

2.5 

3.8 

8.3 

2 

8.1 

3.1 

Readout 

Mode 

drift  scan 

ds/step  stare 

step 

stare 

step  stare 

step  stare 

step  stare 

step  stare 

step  stare 

Exposure 

(sec) 

150 

150 

20 

60 

45 

5 

60 

23 

Revisits 

3 

3 

3 

3 

4 

5 

4 
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2.3  Future  Programs 

As  stated  in  Chapter  1,  two  proposed  future  efforts,  Pan-Starrs  and  LSST,  were 
detailed  in  the  AoA  submitted  to  Congress  by  NASA  in  2007  [3].  Since  the  submittal  of 
the  AoA,  plans  for  both  future  programs  have  matured  significantly.  The  LSST 
scheduled  first  light  is  in  2018  and  no  hardware  has  yet  been  delivered  [10].  The  first  of 
four  identical  pieces  of  Pan-Starrs,  PS  1  has  been  installed  on  Haleakala,  in  Hawaii  and 
has  been  officially  conducting  scientific  observations  since  May  2010.  PS1  will  be  used 
to  test  the  design  and  technology  being  developed  for  Pan-Starrs.  The  second  piece,  PS2, 
is  scheduled  for  installation  in  early  2013.  The  installation  of  the  remaining  pieces  is 
currently  unscheduled  [11].  Pan-Starrs  has  implemented  the  first  major  change  in  NEO 
detection  methods,  and  no  longer  looks  for  moving  objects  based  on  two  single  images 
that  are  both  susceptable  to  noise.  Instead,  Pan-Starrs  creates  a  running  average,  or 
"master  image"  of  the  sky  and  uses  it  in  an  image  differencing  algorithm.  The 
availability  of  the  master  image  is  the  genesis  of  this  research  effort  and  for  this  reason 
the  remainder  of  the  focus  of  this  section  will  be  the  review  of  literature  as  it  pertains  to 
the  methods  and  capabilities  of  Pan-Starrs  and  how  the  Pan-Starrs  program  conducts 
NEO  surveillence. 

2.4  Pan-Starrs  Specifications 

PS1  has  a  1.8  meter  diameter  concave  primary  mirror  and  an  effective  focal 
length  of  8  meters.  The  camera  on  the  PS1  is  an  8  by  8  array  of  CCD  devices 
approximately  5  cm"  each.  The  individual  devices  are  made  up  of  an  8  by  8  array  of 
CCD  cells  with  almost  600  by  600  pixels  per  cell,  providing  approximately  1.4  Giga- 
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Video  output  Analog  out 


Figure  4.  Pan-Starrs  CCD  Structure  [11] 

pixels  total,  as  illustrated  in  Figure  4.  The  pixel  resolution  of  the  10  pm  pixels  is 
approximately  0.258  arcsec  [1 1]. 

PS1  can  utilize  one  of  several  filters  depending  on  the  application.  For  NEO 
observations,  a  wideband  0.5  to  0.8  pm  filter  will  typically  be  used.  With  the  wideband 
filter,  sky  background  noise  is  expected  to  be  approximately  7  electrons  per  pixel,  while 
the  read  noise  in  the  CCD  cells  is  expected  to  be  about  5  electrons.  Typical  exposure 
times  will  be  30  seconds,  in  which  about  seven  square  degrees  of  the  sky  will  be  imaged 
[11]. 

When  the  full  Pan-Starrs  is  operational,  images  from  each  of  the  four  telescopes 
will  be  compared  and  combined  into  a  composite  image  providing  resilience  to  cosmic 
rays,  error  due  to  gaps  between  CCD  cells,  and  bad  pixels.  Any  composite  image  that  is 
found  to  have  no  objects  of  interest  will  be  used  to  build  up  the  running  master  image  of 
the  sky.  New  composite  images  will  be  analyzed  for  objects  above  a  certain  intensity 
threshold,  which  will  be  stored  in  a  database  for  future  reference.  New  images  may 
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optionally  be  smoothed  with  a  modeled  Point  Spread  Function  (PSF)  estimated  from  the 
image  and  will  be  subtracted  from  the  cumulative  master  image  leaving  only  images  of 
objects  that  have  moved  or  changed  intensity  [12].  The  images  remaining  will  be  added 
to  a  separate  database  for  further  review.  Finally,  the  composite  image  is  discarded  to 
free  up  data  storage  [11]. 

2.5  Relevant  Research 

Research  conducted  in  2004  and  published  in  the  October  2005  issue  of  The 
Astronomical  Journal,  “Matched  Filter  Processing  for  Asteroid  Detection”  described 
methods  to  implement  a  matched  filter  algorithm  for  the  specific  application  of  near 
Earth  asteroid  detection  [13],  The  results  of  this  paper  showed  a  40  percent  increase  in 
asteroid  detection  rates  compared  to  the  “step  stare”  and  “drift  scan”  algorithms  [13].  A 
paper  published  in  The  Astronomical  Journal  in  2005,  “Likelihood-Based  Method  for 
Detecting  Faint  Moving  Objects”  develops  a  Maximum  Likelihood  Estimate  of  the 
magnitude  and  velocity  of  dim  moving  objects  given  the  photon  distribution  received 
from  a  5  minute  exposure  time  [14].  The  results  of  this  research  showed  improvement  in 
the  detectable  magnitude  of  objects  for  a  given  telescope,  but  the  process  requires 
extremely  long  exposure  times,  multiple  images,  and  a  significant  amount  of 
computational  time  [14], 

Another  research  effort,  “Detecting  Near-Earth  Objects  Using  Cross-Correlation 
with  a  Point  Spread  Function”  was  conducted  at  AFIT  and  focused  on  the  Linear  program 
telescope  [15].  The  results  of  this  research  found  that  the  Linear  program  would  improve 
detection  by  using  electronic  matched  filtering  and  proper  sampling  [15].  Pan-Starrs  has 
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an  optional  setting  to  smooth  the  difference  image  with  a  modeled  PSF,  in  other  words  to 
perfonn  matched  filtering.  Sampling  turns  out  not  to  be  an  issue  with  Pan-Starrs  as  it 
was  with  Linear.  The  Nyquist  sampling  theorem  dictates  that  the  Pan-Starrs  detector  is 
sampled  as 


A  < 


Xf  ^0.5 x  10  6 meters j ( timet ers ) 


2D 


2(1.8  meters} 


=  1.11  jum , 


(1) 


where  As  is  the  sample  spacing,  A  is  the  wavelength,  /  is  the  focal  length,  and  D  is  the 
aperture  diameter  [16].  Although  it  appears  that  the  10 jum  physical  pixel  size  is  under 
sampling  by  a  factor  of  nearly  10,  the  cutoff  frequency  due  to  the  atmosphere  is  much 
lower.  Fried’s  seeing  parameter,  r0 ,  is  a  measure  of  the  strength  of  the  turbulence  in  the 
atmosphere  and  has  the  effect  of  limiting  the  effective  diameter  of  the  system,  if  the 
diameter  is  larger  than  r0 .  As  shown  in  Figure  5,  the  average  r0  at  Haleakala,  the 
location  of  PS1,  is  approximately  9  cm  during  hours  of  darkness  [17]. 


Figure  5.  Average  hourly  r0  for  Haleakala,  May-Oct 
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This  means  that  r0  is  the  limiting  factor  in  the  spatial  resolution  of  the  PS  1  system,  and 

using  9  cm  for  I)  in  Eq.  (1)  yields  a  minimum  effective  spatial  sampling  of 

;  f  (0.5x\0  6 meters)(8meters) 

As  <  =  - - - - }—x - -  =  22.2 /um  , 

2r0  2(0.09 meters) 

which  means  that  the  10/r/u  physical  pixel  size  of  PS1  is  sampling  at  a  rate  about  twice 
that  which  is  necessary  for  the  average  atmosphere  and  is  fine  enough  sampling  for  an  r0 
of  up  to  20  cm. 

2.6  Summary 

Six  of  the  largest  NEO  detection  programs  over  the  last  25  years  all  use  the  same 
fundamental  approach.  The  newest  program,  Pan-Starrs,  introduces  the  method  of 
producing  an  average  image  of  the  sky  and  using  that  image  for  difference  detection  in 
order  to  identify  new  objects  that  are  not  typically  in  that  image.  Not  much  has  been 
reported  about  the  LSST  program  due  to  the  stage  of  development.  LSST  plans  on  using 
both  the  step-stare  and  image  difference  detection.  Instead  of  building  their  own  master 
image,  the  plan  calls  for  a  star  chart  used  in  a  similar  manner  as  Pan-Starr's  master  image 
[10].  The  availability  of  Pan-Starr's  master  image  provides  a  statistical  tool,  the  expected 
value  of  the  sky  image,  which  may  be  utilized  in  a  classical  signal  detection  approach 
using  a  Likelihood  Ratio  Test,  a  method  that  has  not  been  applied  to  the  NEO  detection 
problem.  The  proposed  algorithm  removes  assumptions  about  the  statistical  nature  of 
photon  arrival  and  implements  a  LRT  detection  scheme.  The  new  algorithm  requires  the 
master  image,  but  uses  it  as  the  statistical  mean  of  an  image  taken  with  the  Pan-Starrs 
camera  as  opposed  to  generating  a  difference  image. 
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3  Methodology 


3.1  Chapter  Overview 

The  results  of  [15]  showed  that  a  matched  filter,  or  a  cross-correlation  with  the 
PSF,  would  improve  NEO  detection  using  the  Linear  telescope.  This  chapter  will  show 
that  using  the  PSF  smoothing  option  of  the  Pan-Starrs  system  on  the  differenced  image  is 
actually  applying  an  assumption  that  the  photon  noise  is  both  white  and  Gaussian. 
Eliminating  this  assumption  is  the  basis  for  the  research  effort.  A  detailed  explanation 
follows  that  describes  the  analytical,  simulation,  and  measurement  tools  used  to  develop 
and  demonstrate  a  more  effective  detection  algorithm  based  on  photons  having  a  Poisson 
distribution.  The  analytical  tools  include  a  derivation  of  the  signal  to  noise  ratio  (SNR) 
as  a  metric  to  compare  the  algorithms.  Simulation  tools  include  the  development  of  the 
complete  model  of  the  Pan-Starrs  equipment,  atmosphere  and  astronomical  objects. 
Measured  tools  include  the  development  of  a  relative  test  environment  to  recreate 
conditions  that  will  be  seen  in  practical  applications  of  the  algorithms. 

3.2  Photon  Statistics 

Detecting  the  presence  of  an  object  optically  is  done  by  measuring  the  number  of 
photons  received  at  the  detector  from  the  object.  There  will  also  be  photons  arriving 
from  the  ambient  background  radiation  of  the  scene.  Detecting  the  object  requires  a 
statistical  characterization  of  the  arrival  of  all  photons  in  order  to  distinguish  between 
photons  from  a  target  object  and  those  from  background  radiation.  Assume  a  total  of  N 
photons  are  emitted  over  time  in  any  one  direction  from  an  isotropic  source.  The  rate, 
A  (t),  at  which  n  photons  arrive  in  t  seconds  at  the  detector  would  be 
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A(t)=,±(p,wtons/s ). 

The  probability,  p  of  a  photon  arriving  in  a  time  interval  (h,t2)  is  simply  [18] 

P=  PaiM*. 

•'*1 

where  pA  (f)  is  the  probability  density,  or  the  photon  arrival  probability  per  second. 
Intuitively,  the  probability  density  is  proportional  to  the  arrival  rate  for  photons,  or 

Pa  {t)ccA{t). 

The  probability  density  must  integrate  to  one;  therefore, 

P  [photon  arriving] 


Pa  (0  =  — 
U  J  Nt 


or 


Pa  (r)  = 


A (t)  f  P [photon  arriving] 


N 

V  / 

The  arrival  of  KN  photons  from  a  source  of  N  photons  in  a  time  interval  (t,  ,t2 ) 
can  be  characterized  as  a  binomial  random  variable  with  N  total  experiments,  KN 
successes,  and  a  probability  of  success  of  p.  The  probability,  P(Kn  )  of  KN  photons 
arriving  in  the  interval  {tl,t2)  is  then, 


P(K„)  =  \N]pk"(1-pTK" 


N\ 


(2) 
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N (N  -i)- ...-(N  - Kn  +l) 


Kn\N> 


J  A{t^dt  1--^- J ' A(t}d\ 


The  total  number  of  photons  emitted  from  a  source,  N  can  be  assumed  to  be  uncountable 
and  large;  therefore, 


°[  kn\n 


(3) 


and  when  N  is  large,  [19] 


>(')<* 


'-ir*')* 


-J  2  A(t)dt 


Eq.  (3)  simplifies  to 


N(N-1)-...-(N-Kn  +l) 


J  A{t^dt 


-ik„  _r- 


l24')dt 


kn\ 


| 1 A  (t)dt 


-  f  2  A (t)dt 

Jtt 


(4) 


The  average  number  of  photons,  KN  arriving  from  a  source  of  N  photons  in  the  time 
interval  (t{,t2)  is 


KN=^  A{t)df, 


(5) 


therefore,  Eq.  (4)  becomes 


P(KU) 


V  KN  p~kn 
^ N  c 

Kn\ 


(6) 
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which  is  the  Probability  Mass  Function  for  a  Poisson  random  variable.  This  random 
variable  is  characterized  by  a  single  parameter  which  is  equal  to  both  the  mean  and  the 
variance  [20], 

The  range  between  the  detector  and  the  object,  and  the  source  of  the  photons  from 
an  object  detennines  how  the  expected  number  of  photons  arriving  at  the  detector  is 
estimated.  A  NEO,  within  1.3  AUs  of  Earth,  is  relatively  close  compared  to  a  distant  star 
light-years  away.  NEOs  are  also  not  the  source  of  photons;  rather  they  reflect  photons 
from  other  sources.  The  expected  number  of  photons,  KN  ,  received  and  converted  into 
photo-electrons  by  the  detector  from  a  NEO  illuminated  only  by  natural  light  can  be 
calculated  using  Eq.  (7)  [16]. 


ts  S IBAxABptijz 'at pD  A  At 
N  4  R2hc 


(7) 


where: 


•  SIB  is  the  solar  spectral  irradiance  incident  on  the  NEO  in  Watts  (J/s)  per  square  meter 
per  micro-meter, 

•  Ax  is  the  bandwidth  of  the  filter  in  micro-meters, 

•  Ab  is  the  two-dimensional  surface  area  of  the  NEO  that  is  nonnal  to  the  optical  axis  in 
square  meters, 

•  pt  is  the  dimensionless  reflectivity  of  the  NEO, 

•  77  is  the  quantum  efficiency  of  the  detector  in  photo-electrons  per  photon, 

•  Ta  is  the  dimensionless  transmittance  of  the  atmosphere, 

•  T0  is  the  dimensionless  transmittance  of  the  optical  system, 

•  D  is  the  diameter  of  the  receiver  aperture  in  meters, 
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•  X  is  the  photon  wavelength  in  meters, 

•  At  is  the  exposure  time  in  seconds, 

•  R  is  the  range  to  the  NEO  in  meters, 

•  h  is  Planck's  constant  in  Joule  seconds  per  photon,  and 

•  c  is  the  speed  of  light  in  meters  per  second. 

The  natural  light  incident  on  the  NEO  is  almost  entirely  from  the  Sun.  The  solar 

spectral  irradiance  at  the  top  of  Earth's  atmosphere,  shown  in  Figure  6,  is  used  because  it 
closely  predicts  the  irradiance  incident  on  a  NEO  without  an  atmosphere  in  an  orbit  close 
to  the  Earth’s.  The  bandwidth  of  the  wideband  filter  typically  used  for  NEO  detection  by 
Pan-Starrs  is  0.3  pm;  however,  this  model  will  simulate  the  monochromatic  response 
over  the  range  of  the  filter  in  0.1  pm  increments,  so  0.1  pm  is  used  for  Ax,  and 


Wavelength  (nm) 


Figure  6.  Solar  Radiation  Spectrum  [21] 

0.5  pm,  0.6  pm,  0.7  pm,  and  0.8  pm  will  be  used  for  X .  Four  wavelength  dependent 
values  corresponding  to  each  X  found  in  Figure  6  will  be  used  for  SIB .  The  irregular  and 
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unknown  shape  of  a  NEO  makes  using  exact  values  for  AB  impossible;  therefore  a 
circular  surface  will  be  assumed  and  the  equation  for  the  area  of  a  circle  with  a  diameter 
greater  than  or  equal  to  140  meters  will  determine  the  value  used  for  AB .  This  is  a 

reasonable  assumption  because  the  variation  in  angular  size  on  the  detector  plane 
between  an  irregularly  shaped  object  and  a  spherical  object  is  not  measurable.  NEOs  are 
expected  to  have  a  reflectivity  greater  than  3  percent,  which  provides  a  lower  bound  on 
the  value  of  p ,  [2].  A  study  done  by  the  Department  of  Physics  at  Harvard  University 

measured  the  total  effective  throughput  of  Pan-Starrs  using  a  calibrated  Silicon 
photodiode  and  a  tunable  laser.  At  even  intervals  of  the  wavelengths  of  the  wideband 
filter,  0.5,  0.6,  0.7  and  0.8  pm,  the  study  found  that  the  optical  transmittances  for  PS1  are 
approximately  0.67,  0.85,  0.99,  and  0.95,  respectively  [22].  These  values  will  be  used  as 
the  values  for  the  product  rj t0  .  The  typical  exposure  time,  At ,  of  a  Pan-Starrs  image 

will  be  about  30  seconds  [11].  The  furthest  range,  R  ?  defining  a  NEO  is  1.3  AUs,  or 
approximately  1.9448  xl0nm  [1].  A  complete  derivation  of  Eq.  (7)  is  found  in  [16]. 

The  expected  number  of  photons  incident  on  the  aperture  over  the  exposure  time 
of  the  detector  for  objects  making  up  the  background  radiation  requires  a  different 
method  than  Eq.  (7).  The  AB  apparent  magnitude  system  is  a  measure  of  the  relative 
irradiance  of  celestial  objects.  Alpha  Lyr  (Vega)  is  the  reference  magnitude  (zero 
magnitude)  on  the  AB  system.  The  number  of  photons  received  from  the  object  is 
proportional  to  the  irradiance  of  the  object;  therefore,  if  the  expected  number  of  photons 
from  Vega  is  calculated,  the  AB  system  can  be  used  as  a  measure  of  the  relative  number 
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of  photons  between  two  objects  [23].  In  the  AB  system,  the  magnitude  of  an  object,  Ms 
is 

M,  =-2.5  log10(F)-48.585,  (8) 


eroS 

where  F  is  the  flux,  in  - - .  The  relative  difference  in  magnitude  between  two  objects, 

cm~ 

Ms  —Mn  =Mx  equates  to  a  factor  of  2.5 12M’ decrease  in  irradiance  or  photons.  Using 

this,  the  relation  between  the  number  of  expected  photons  received  between  an  object  of 
interest  and  the  reference  Vega  is 

M  —Mv  —M  , 

s  V  s' 


2.512Mi  =i, 
K. 


K=  ^ 


s  2.5 12"*  ' 


(9) 


where  Mv  =  0,  the  magnitude  of  Vega,  Kv  is  the  expected  number  of  photons  received 
from  Vega  at  the  detector,  Ks  is  the  expected  number  of  photons  received  from  the 
source,  and  Ms  is  the  magnitude  of  the  source,  say  a  star,  on  the  AB  magnitude  system 
and  can  typically  be  found  from  a  look-up  table.  Solving  Eq.  (9)  requires  Kv  ,  found  by 


Ky  = 


F  B  A  - At 


(10) 


where: 

•  F  is  the  flux  in  Joules  per  square  meter, 

•  B  is  the  bandwidth  of  the  Pan-Starrs  wideband  filter  used  in  NEO  detection  in  Hz, 
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•  A  is  the  area  of  the  aperture  in  square  meters, 

•  At  is  the  exposure  time  in  seconds,  and 

•  E  is  the  energy  per  photon  in  Joules  per  photon. 

To  find  the  flux,  set  M s  =  Mv  =  0  and  solve  for  F  in  Eq.  (8),  and  then  convert  the  units 


of  the  flux  to  equal  the  units  of  the  flux  in  Eq.  (10): 

-20  ergs 


F  = 


3.6813x10" 
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10  !  J 

V  erS  j 


1  a4  2 

10  cm 
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=  3.6813 xlO"23  -4. 
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(11) 


The  energy  per  photon  is  wavelength  dependent,  and  the  wideband  filter  passes  500  mn 
to  800  mn  wavelengths.  A  close  approximation  of  the  expected  number  of  photons 
received  will  sum  the  number  of  photons  received  for  each  of  the  wavelengths  500  mn, 
600  mn,  700  mn,  and  800  nm.  The  calculation  for  each  of  the  wavelengths  will  divide 
the  bandwidth  of  the  wideband  filter  into  four  equal  bands:  500-575  nm,  575-650  nm, 
650-725  nm,  and  725-800  nm.  At  500  nm,  each  photon  has  a  frequency,  /  of  roughly 

6.0xl014/fe.  The  energy  in  each  photon  at  500  nm,  Z?,is 

E  =  hf 

=  (6.626068x10"34/-s)(6.0x1014Hz) 

=  3.9756x10 ~19/ 

where  h  is  Plank's  constant.  The  bandwidth,  B  ,  used  in  Eq.  (10)  for  this  wavelength  is: 

B  =  — - - - — 

500nm  51 5  mn 

=  7.82609  xlO13  Hz 

n 
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where  c  is  the  speed  of  light.  The  number  of  expected  photons  from  Vega  in  the  band 
500  to  575  mn,  Kvl  received  by  the  1.8m  Pan-Starrs  mirror  in  a  30  second  exposure, 


using  Eq.  (10)  is: 


Kv  i=- 


3.6813x10~23^4  (7.82609xl013/7z)(305)(0.9m)2^ 

777  )  '  ' 


3.9756xl(T19 1 


photon 


=  5.53215x10“  photons 

Similar  calculations  for  the  remaining  three  wavelengths  yields  the  number  of  expected 
photons  from  Vega  in  the  band  575  to  650  mn,  Kv ,  =  5. 1 066  x  1  ()' 1  photons  ,  in  the  band 

650  to  725  mn,  Kv3  =4. 72506 x\0n  photons ,  and  in  the  band  725  to  800  mn, 
KV4  -  4.38756x10“  photons .  The  total  number  of  expected  photons  received  in  a  30 

second  exposure  using  the  wideband  filter  of  Pan-Starrs  is  the  sum  of  the  number  of 
photons  from  each  wavelength: 


Ky  —  Kvl  +  Kv  2  +  KV3  +  KV4 , 
=  1 .975 13  x  1012  photons . 


(12) 


Now,  given  the  magnitude  of  a  celestial  object,  Eq.  (9),  and  Eq.  (12),  it  is  possible  to 
calculate  the  expected  number  of  photons  received  by  the  detector  of  Pan-Starrs  from  any 
object.  This  will  be  used  to  model  stars  during  simulations  with  varying  brightness 
which  contribute  varying  amounts  of  background  radiation  in  the  scene  where  a  very  dim 
NEO  may  be  located. 

Using  Kn  from  Eq.  (7),  Kv  from  Eq.  (12),  and  Eq.  (9),  it  is  possible  to  calculate 
the  apparent  magnitude  of  the  NEO.  By  setting  K s  =  K N  and  Ms  =  M N  in  Eq.  (9)  and 
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solving  for  M N  by  taking  the  natural  logarithm  of  both  sides,  the  apparent  magnitude  of 
the  NEO,  Mn  is 


M„  = 


ln(gy)-ln(K„) 

ln(2.512) 


3.3  Likelihood  Ratio  Test 


A  LRT  is  a  detection  test  that  compares  the  ratio  of  two  probabilities  conditioned 
alternately  on  two  hypotheses  to  a  threshold.  A  complete  derivation  of  the  LRT  can  be 
found  in  [16].  The  LRT  is  defined  as 


P(d(x,y)V(x,y)s[l,Md]\Hl)  > ^ 
P(d(x,y)y(x,y)e[l,Md]\H0)<  1 


(14) 


Recall  that  M d  is  the  number  of  pixels  in  one  dimension  of  the  simulated  square  detector 
plane.  Hx  is  the  hypothesis  that  a  NEO  is  present,  and  H0  is  the  hypothesis  that  a  NEO 
is  not  present.  p(yd(x,y)V(x,y)e[l,Mc,]\Hi )  is  the  probability  of  the  data  at  the  (x,y)  pixel 

given  hypothesis  i  e  {0,1} .  The  bounds  on  x  and  y  will  be  assumed  for  the 

remainder  of  the  document  and  the  notation  will  be  shortened  to  P[d(x,  y)  |  H~^ .  In 

practical  applications,  a  Neyman-Pearson  detection  test  in  which  the  detection  threshold 
is  set  to  a  value  that  produces  a  desired  false  alarm  rate  should  be  used  in  NEO  detection. 
For  purposes  of  comparing  the  performance  of  the  current  and  proposed  methods,  a 
Bayesian  detection  test  is  used  in  order  to  sweep  through  all  possible  values  of  thresholds 
and  generate  a  ROC  curve. 
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As  previously  stated,  the  number  of  photons  reflected  from  a  NEO  that  are 
converted  to  photo-electrons  by  the  detector  of  PS  1  has  a  Poisson  distribution  with  mean 
and  variance  equal  to  KN  from  Eq.  (7),  similarly  photons  from  background  radiation  are 


Poisson  with  mean  and  variance  equal  to  Ks  from  Eq.  (9).  Given  the  properties  of 

Poisson  random  variables,  the  composite  image,  or  the  sum  of  the  photons  from  all 
objects  in  the  scene  is  also  Poisson.  The  current  method  of  detection  by  Pan-Starrs 
without  making  simplifying  approximations  about  the  distributions  of  the  photons  does 
not  lend  itself  to  a  LRT.  Subtracting  the  master  image,  which  is  the  deterministic  mean 
of  the  sky  image,  from  one  exposure,  which  is  Poisson  distributed,  results  in  a  random 
variable  whose  mean  and  variance  can  be  calculated,  but  whose  mass  function  is  no 
longer  Poisson.  One  property  of  the  Poisson  random  variable  that  simplifies  this  problem 
is  that  the  Poisson  random  variable  is  well  approximated  as  a  Gaussian  random  variable 
with  mean  and  variance  equal  to  the  parameter  K ,  under  high  illumination  conditions. 
The  difference  between  a  Gaussian  random  variable  and  a  detenninistic  variable  is  still 
Gaussian,  with  a  shifted  mean.  This  approximation  makes  a  LRT  much  more  feasible 
and  its  derivation  follows  [20] . 


/>[</(*, y)|ifi]> 
G  P[j(x,y)|tf0]< 

Ho 


(15) 


where: 

d  |  Hl  ~  Gaus  |  Kb  +  ( KN  +  Ks  )  *  g,  <j2  j 
d  |  H0  ~  Gaus  {kb  +  Ks  *  h,  cr2  j 

•  d  is  the  data, 

•  Hx  is  the  hypothesis  that  a  NEO  is  present  in  the  scene, 
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•  H0  is  the  hypothesis  that  a  NEO  is  not  present  in  the  scene, 

•  Kb  is  the  expected  number  of  photons  received  from  background  noise, 

•  Ks  is  the  expected  number  of  photons  received  from  known  objects  in  the  image, 

•  Kn  is  the  expected  number  of  photons  received  from  a  NEO, 

•  g  is  the  impulse  response  of  the  system, 

•  cr 2  is  the  variance  (set  equal  for  both  hypotheses  to  reproduce  Pan-Starrs  detection 
scheme),  and 

•  77  is  the  threshold  of  detection. 


The  master  image,  I M  is 


I  m  Kb  +  Ks  *  g  , 


(16) 


and 


r  ,  x  -1  t— r-i—r  1  I-—1 ^-[<<(*>)’)_*i>-(Arw+*rs)*sT 

x  y  y  G 


Using  properties  of  convolution, 


x  y  v  CJ 


1  t-v.v)-^**]2} 


(17) 


The  objects  from  which  photons  are  arriving  can  be  viewed  as  a  point  source  (a  scaled 
Dirac  delta  function)  due  to  the  fact  that  the  angular  size  of  the  objects,  ignoring  larger 
objects  within  the  solar  system,  are  smaller  than  one  pixel  in  the  detector  plane;  therefore, 
in  general 

K  *  8  =  s(am  ~x-a,pm-y-  p)g  (a,  0) , 

a  p 

=  Kmg(am-x,pm-y),  (18) 
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where  is  the  location  of  the  point  source  with  magnitude  Km,  the  average 

number  of  photons  received  from  object  m  .  For  the  remainder  of  the  document,  photon 
source  objects  will  be  from  stars  located  at  (as,/?s)  with  an  average  number  of  received 

photons  of  Ks,  NEOs  located  at  ( aN,/3N )  with  an  average  number  of  received  photons 
of  Kn  ,  and  average  background  photons  at  each  pixel,  KB  . 

Now,  Eq.  (16)  becomes 

1M  (*,  y)  =  Kb  +  Ksg  (as  - X,  ps  -  y) ,  (19) 


and  Eq.  (17)  becomes 

P\d(x,  y)l^1]  =  llH  p—- 

X  y  V27TCr 


1(7 


(20) 


similarly, 


p[rf(jp,y)iH0]=nn 


2cr' 


[d(x,y)-IM(x,y)] 


■sJTjtct 


(21) 


Computing  the  LRT  using  Eq.  (20)  and  Eq.  (21),  Eq.  (15)  becomes 

F[rfq,y)|g,]> 

0  P[d{x,y)\H„\< 

Ho 


(22) 


It  can  be  shown  that: 
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n/w 


f{xi)f{x2)---f{-x„)  _f{-x i)  /(*2)  f  M  _yrf  (xi) 
g(yi)g(y2)---g(yn)  g(y i)  sM  sW  /=f  ^ ( y, )  ’ 


Eq.  (22)  simplifies  to 


Ac=nrr 


--^Tp (x,y)-Iu  (. x,y)-KNg(aN  -x,fiN  -y)]'  }  Hi 


^j[dfi,y)-/Mfi,y)]2| 


\--^2[d(x,y)-IM(x,y)-Kflg(aN-x,/3N-y)]2+-yj[d(x,yyiM{x,y)j\  > 

AG=llfle  J  rj . 

x  y  ^ 


Next,  square  the  terms  of  the  exponential,  combine  like  terms  and  simplify: 


A0=nn< 


~y\  --Z-8  iaN -x>Pn -y)+d(x^y)g(^N ~x’Pn  -y)-lM  ( x,y)g(aN  -x,pN -y)  y  > 


Ac  =e 


Y,Yj^Y^g2^a«^x^N^y)+d(x’y)s(aN-x,/}N-y)-IM(x,y)g{aN-x,/}N-y)^ 


Finally,  take  the  Natural  Logarithm  of  both  sides  and  simplify,  noting  that  — >  o : 


hiAG='E'Ld(x’y)s(aN-x’PN-y)-'E'LIM{x’y)g(aN-x’PN!My)>  ^+IIy«2(a»-^A-3')’ 

*  y  x  y  5  ^ n  x  y  ^ 

H 0 

where  the  right  side  of  the  equation  is  some  constant,  y  and  ( aN ,  f3N )  are  the 
hypothesized  coordinates  of  the  NEO,  and  fall  within  some  range  («,/?)  e  [l,  ] . 


biAG(a,/?)  =  ^]^tf(x,y)g(a-x,^-y)-^^/M(x,y)g(a-x,/?-y)  y, 

XV  XV 


27 


=  {d~IM)*g><Y  ■ 

«o 


(23) 


Eq.  (23),  which  precisely  describes  the  Pan-Starrs  method,  shows  that  using  peak 
detection  on  the  convolution  of  the  difference  image  with  the  PSF  as  a  detection 
algorithm,  has  built  into  it  the  assumption  that  the  photon  data  received  is  Gaussian  and 
that  the  variance  of  the  photon  data  is  equal  for  both  H ]  and  H0 . 

If  a  Poisson  model  is  used  for  the  noise  in  the  imagery,  the  log-likelihood  ratio 

test  is: 


Ap  = 


P[d(x,y)\Hl~\  > 
P[d(x,y)\H0~\  < 


n , 


(24) 


where, 


p[X*,;>ytf,]=nn 


Using  Eq.  ( 1 6)  and  Eq.  (18): 


d\Hi~Pois{KB+(KN  +  Ks)*g} 
d\H0~  Pois{RB  +  Ks*g} 

[KB+(Kx+Ksygf-\ 


d  (v,  y)! 


[i  (x  y)  +  K  s(a  -xB  _  {~^B,y)-KNg(aN-xjN-y)} 

p[rf(x.y)iH1]=nnL'“  'J. — - .  (25) 


x  y 


d{x,y)\ 


similarly, 
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(26) 


'Ww)iB.]=nn 


x  y 


d(x,y)\ 


Computing  the  LRT  using  Eq.  (25)  and  Eq.  (26),  Eq.  (24)  becomes 


Jm  (x’  y)  +  KNs{aN-x,pN-y)_ 


d(x’y )  A-Iu(x,y)-K!ig(alf-x,filf-y)} 


A,=nn 


dpzryfi 


d^  e{-iM(*.y)} 


H, 

> 

< 


V 


dj^rfy. 


Ac 


nn 


IM(x,y)  +  KNg(aN-x,^N-y) 
!M{x,y ) 


d(x<y)  I  -KNY^s(aN-x,pN-y)-Y^JjA<y) 

e 


Take  the  natural  logarithm  of  both  sides  and  simplify, 

<  ^ N  x  y 

Ho 

The  threshold,  ij ,  for  Bayesian  detection  where  the  prior  probability  of  a  NEO  being 
present  is  equal  and  unifonn  costs  are  assumed,  is  set  equal  to  one.  It  is  assumed  that  the 
prior  probability  of  a  NEO  being  present  or  not  is  equal,  which  would  be  grossly 
inaccurate;  however,  the  goal  is  to  compare  the  perfonnance  with  the  use  of  ROC  curves, 
which  vary  the  threshold  over  all  possible  values.  The  number  of  photons  received  from 
a  NEO,  Kn  is  unknown.  The  simplifying  assumption  on  the  prior  probability  of  the 

presence  of  a  NEO,  allows  the  tenn  on  the  right  side  of  the  equation  containing  KN  to  go 
to  zero.  It  will  be  shown  in  the  next  section  that  KN  ,  when  estimated  depends  on  the 
data  received,  and  therefore  should  not  be  on  the  right  side  of  the  equation.  Without  any 


d  (x,  y)  |n  L  +  KNg  {aN-x,PN-y)\ 


iu(*>y) 
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usable  data  to  formulate  a  value  for  // ,  the  inaccurate  assumption  to  set  it  to  1  is 
necessary.  The  remainder  of  the  right  side  of  the  equation  will  be  set  to  some  constant, 
y  in  order  to  sweep  through  all  possible  thresholds  and  create  a  ROC  curve.  As  with  the 

Gaussian  LRT  derivation,  ( aN,/3N )  are  the  hypothesized  coordinates  of  the  NEO,  and 
falls  within  some  range,  (a,/?)  e  [l,Mrf] . 


lnA,(a,/?)  =  XZ^k>- 


KNg(a-x,p-y )  > 

1+ - n — t - r  r 


(27) 


It  should  be  noted  that  Eq.  (27)  does  not  estimate  the  location  of  a  NEO,  it  merely  detects 
the  presence  of  a  NEO  in  the  image.  Estimation  theory  can  be  applied  to  estimate 
(aN,/3N),  the  simplest,  but  not  necessarily  the  optimal  estimate  would  be  a  peak 

detection  on  the  two  dimensional  likelihood  matrix,  In  Ap  ( a,j3 ) . 


3.4  NEO  Magnitude  Estimation 

The  total  number  of  photons  received  in  an  image,  K  is  the  summation  of  the 
number  of  photons  received  at  each  pixel,  d 

k =XX^(*,>0-  (2g) 

*  y 

The  summation  of  Poisson  random  variables  is  Poisson  distributed  with  mean  and 
variance  equal  to  the  sum  of  the  means  and  variances  of  each  component  random  variable 
[20] .  The  mean  of  K  is  the  summation  of  the  expected  number  of  photons  received  by 
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each  source,  in  this  case  the  photons  received  by  the  NEO,  KN  ,  the  photons  received  by 
any  background  stars,  Ks ,  and  the  photons  received  by  background  radiation,  KB 

K  =  Kn+Ks+KbM],  (29) 

where  M d  is  the  number  of  pixels  in  one  dimension  of  the  square  detector  plane.  Given 

that  the  PSF  is  normalized,  performing  a  summation  of  all  of  the  pixels  of  the  master 
image  using  Eq.  (19)  yields, 

YJHIM(x’y)  =  YjH[^B+Ksg(as~x,/3s-y)\, 

x  y  x  y 

=  KbM]+Ks.  (30) 

A  Generalized  LRT  (GLRT)  is  formed  when  an  unknown  parameter  is  replaced  by  the 
maximum  likelihood  estimate  for  that  parameter  [24].  The  maximum  likelihood  estimate 

of  Kn  ,  Kn  can  be  found  by 

Kn  =  argmax{lnP(if  |  ^)}  , 


:  arg  max  <; 

kn 


4+XS/«(h.y) 


In 


x  y 


Kl 


=  arg  max  <K  In 

Kn  i 


Kn  (x’-y)  ~Kn  “ZZ7m  ■ 
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The  argmax  can  be  found  by  setting  the  partial  derivative  of  the  log-likelihood  with 
respect  to  KN  equal  to  zero, 

o=A- rta  *»+ZE/«tt-y)  • 

[  v  x  ?  y  *  y  J 

=  ^ _ x _ i, 

x  y 

KN=K-YjYjIm{X^)' 

x  y 

=ZZ^/(x’>;)_ZE/M(hi)-  (3i) 

*  y  i  j 

Using  this  estimate  in  the  implementation  of  the  Poisson  GLRT  algorithm  not  only 
introduces  additional  randomness  into  the  GLRT,  but  it  is  dependent  on  the  master 
image.  The  effect  of  this  estimate  will  vary  with  the  magnitude  of  any  background  stars 
in  the  image.  The  mean  of  the  estimate  is 

E  kN  =Kn  (32) 

The  variance  of  the  estimate  is 

Var(kN)  =  KN+^^IM(x,y)  (33) 

'  /  x  y 

This  shows  that  the  estimator  is  unbiased  and  the  variance  of  the  estimate  will  increase 
with  the  increase  in  photon  counts  from  the  master  image.  Analytically,  the  signal  to 
noise  ratio  (SNR)  will  be  used  to  predict  whether  the  Poisson-based  GLRT  or  the 
Gaussian-based  LRT  will  perform  best.  The  SNR  will  be  defined  as 
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SNR 


(34) 


E 

1  1 

a 

J; 

H\ 

-E 

A  (  Ct-N  1  Pn  ) 

< 

>/ 

War(N(aN,PN)) 

This  definition  accounts  for  the  magnitude  of  separation  of  the  mean  between  the  two 
hypotheses.  This  reveals  the  ability  of  the  algorithm  to  distinguish  the  presence  of  a  NEO 
as  the  LRT  is  compared  to  a  threshold  while  also  not  making  false  positive  detections. 
The  further  the  value  of  the  LRT  given  H l  is  from  the  value  of  the  LRT  given  H0 ,  the 

more  detectable  the  object  is  without  false  positives.  The  statistical  fluctuation  of  the 
LRT  is  indicative  of  the  ability  of  the  noise  to  drive  the  LRT  value  above  or  below  the 
threshold.  This  definition  of  SNR  assumes  symmetry  conditions  with  the  variance  of  the 
LRT.  While  the  variance  of  the  LRT  given  the  two  hypotheses  is  not  equal,  the  variance 
given  Hl  is  treated  as  an  upper  bound  on  the  unconditioned  LRT  and  the  variance  of  the 
LRT  given  H0  is  treated  as  a  lower  bound  on  the  unconditioned  LRT  enabling  the  use  of 

Eq.  (34)  as  the  definition  of  SNR.  The  SNR  of  the  LRT  as  defined  is  a  metric  that 
reveals  the  ability  of  the  algorithm  to  achieve  a  high  probability  of  detection  and  a  low 
probability  of  false  alarm. 

3.5  SNR  Analysis 

The  following  analysis  compares  whether  the  SNR  of  the  Poisson  GLRT  is  higher 
than  the  Gaussian-based  LRT,  the  larger  SNR  will  be  indicative  of  the  better  perfonning 
algorithm.  First,  the  elements  of  Eq.  (34)  are  derived  for  the  Poisson-based  GLRT. 
Looking  at  the  logarithmic  element  of  Eq.  (27),  the  first  approximation  of  the  Taylor 
Series  expansion  is 
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(35) 


i  f,  KNs{ccN-x,pN-y)\  KNg(aN-x,pN-y) 

I  1  M*>)0 

for  _1<g„g(g  -x,A-y)<1 

Mwv) 

A  photon  count  of  less  than  zero  does  not  make  physical  sense;  therefore,  this 
approximation  holds  when 


r?  M*»y) 

g(aN-x,pN-y) 

Using  the  simulated  PSF  for  Pan-Starrs  and  ranging  the  magnitudes  for  the  background 
stars  from  the  dimmest  measureable  star  to  the  brightest  star,  the  required  value  for  K N 
for  this  approximation  to  hold,  regardless  of  the  magnitude  of  the  background  image  was 
found  to  be 


<15, 000  photons  . 

Eq.  (33)  shows  that  the  variance  of  the  maximum  likelihood  estimate  for  KN  may  cause 

this  restriction  to  be  exceeded  when  the  background  image  contains  a  large  number  of 
photons.  When  this  scenario  is  encountered,  the  following  SNR  analysis  no  longer 
accurately  predicts  the  perfonnance  of  the  LRT. 


a  y  K., 


t  |  KNg(a-x,p-y) 


d{x,y)KNg{a-x,p-y) 

KNIM(x,y) 


The  expected  value  of  the  Poisson-based  LRT  given  Hx  is: 

,d(x,y)g(a-x,/3-y) 


Ap(a,j3)  Hr 


4  (u  >’) 
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=zz 

x  y 


E[cl(x,y)\g{a-  x,  p-y ) 

Mat) 

{KNg(aN  - x’ Pn  - y) g(a - X, p - y)  + 1 M  {x,y)g(a-x,p-y)) 


rM(x,y) 


■i+ZZ 


KN82{aN-x>PN-y ) 


*  y 


The  expected  value  of  the  Poisson-based  LRT  given  H0  is: 

E  [_d  (x,  y  )]  g(a-x,p~  y) 


(36) 


£[A,(a,/?)|H0]«XZ: 


x  y 


JM  {x,  y ) 


^^jjJc<y)g{a-x,p-y)  _1 


V(<y 


The  variance  of  the  Poisson-based  LRT  given  is: 


(37) 


Var(^Kp^a,P)\Hxy. 


V  *  T  V'v?  J  )  J 
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ZZ 


d(x,y)g(a-x,P~y) 


'Till 

xi  y\  xi  yi 


JM{x,y) 

E[d(xl,yl)~]E[d(x2,y2)']g(a-xl,p-yl)g(a-x2,p-y2) 

I M  (Xl>  Jl  )  (-L  ’  ^2  ) 

(£'[<i(x1,y1)])  ^(a-Xpjff-y,) 

(A’  fl  ) 

j£'[rf(x,y)]  +  (£'[t/(xpy1)])2Jg2(a-x,,/?-y) 


-zz 
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KNg2{aN-x,PN-y) 

JM(x,y) 
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*1  y\  x2  y2 


"£[d(x1,y1)]§(«-x1,/?- Vj)^ 

^£[d(x2,y2)]g(a-x2,/?-y2)^ 

,  M^i)  , 

v  I M  (x2,  y2  )  j 
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( E[d(x,y)\g2(a-x,p-y ))  ^  [E[d (x„ yj 

l2M(x,y) 
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Ksh(x- 


aN-x,PN-y) 
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V 
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([^ygKr 


-x,PN-y)g2(a-x, 
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i  m  (x>y) 


Var(AP(aN,pN)\Hl)*YJ'E, 


KNg3(aN-x,pN-y) 
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The  variance  of  the  Poisson-based  GLRT  given  //0  is: 


(38) 


Var(KP{a,p)\Ha)^YjYjYjYj 


X\  Ti  x2  y2 


E[d(xl,yl)\g{a-xl,p-yM(  E[d(x2,y2)]g(a-x2,p-y2) 


MWi) 
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hA^yP) 
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Iu{xx,yx)g(a-xl,p-y1)\(  L 


I m  {x2, y2) 
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IM{x,y)g2{a-x,P-y) 
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Var  ( Ap  ( aN ,  pN  )|  H0 )  *  £  £ 


g2(aN-x,PN-y ) 

Lix^y) 
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The  difference  between  Eq.  (38)  and  Eq.  (39)  is  the  tenn 


EE 


KNg\aN-x,PN-y) 


x  y 


which  is  much  less  than  Eq.  (39)  for  small  diameter  NEOs.  The  SNR  as  defined  by  Eq. 
(34)  is  bounded  by 


zz 


Nvtr  («,v  -x,pN-y) 

G(x,}’) 


zz 


(a;v  x>Pn  x>Pn  y) 
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g2(aN-x,PN-y) 
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where  the  bounds  are  approximately  equal  for  small  NEOs.  Eq.  (40)  is  independent  of 


the  maximum  likelihood  estimate,  K N  ,  which  suggests  that  the  value  used  for  KN  does 


not  affect  the  performance  of  the  LRT.  Recall  that  Eq.  (40)  was  derived  with  the 
approximation  in  Eq.  (35),  and  that  the  conditions  under  which  this  approximation  holds 

can  be  exceeded  if  the  value  used  for  KN  is  large.  When  the  variance  of  this  estimate 


causes  the  restrictions  on  the  approximation  in  Eq.  (35)  to  be  increasingly  exceeded,  the 
net  effect  of  using  this  estimate  is  an  increased  variance  of  the  value  of  the  GLRT.  The 
increased  variance  of  the  GLRT  is  due  to  one  of  the  uses  of  the  estimate  being  inside  a 
logarithmic  tenn  and  the  other  being  in  the  denominator  outside  the  logarithmic  tenn. 
The  increased  variance  of  the  GLRT  will  cause  a  decrease  in  the  expected  SNR.  Given 

that  the  SNR  does  not  depend  on  the  value  for  KN  unless  it  is  large,  and  the  variance  of 


Kn  can  be  very  large  depending  on  the  background  master  image,  an  easy  way  to 

eliminate  the  possibility  of  KN  being  inaccurately  large  is  to  set  it  equal  to  a  small 
constant  rather  than  its  maximum  likelihood  estimate.  Eq.  (40)  predicts  this  approach 
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will  yield  the  best  performance  for  the  Poisson  LRT  under  all  environmental  scenarios; 
therefore,  the  Poisson-based  LRT  as  opposed  to  the  GLRT  will  be  used  with  the 

unknown  parameter  KN  set  equal  to  one.  Simulations  will  be  used  to  verify  this 
qualitative  analysis. 

A  similar  derivation  produces  the  SNR  for  the  Gaussian-based  Pan-Starrs  method. 
The  expected  value  of  the  Gaussian-based  LRT  given  Hx  is: 

E[AG(a^)\Hi~\  =  YjYj{ELd(x’  30]#  (a-x,fi-y)-IM(x,y)g(a-x,f3-y)}, 

x  y 

KNg(aN-x,pN-y)g(cc-x,p-y) 

x  y 

E  [  Ac  (aN,PN)\  J  ]  =  XX  Ksg  2{aN-x,PN-y).  (41) 

x  y 

The  expected  value  of  the  Gaussian-based  LRT  given  H0  is: 

E[_AG(a,P)\Ho  ]  =  YjYj{E[d(x’y)]s{^-x,p-y)-IM(x,y)g(a-x,p-y)}, 


x  y 


£[AG(a,/?)|tf0]  =  0.  (42) 

The  variance  of  the  Gaussian-based  LRT  given  Hx  is: 

Var ( Ag  (a,p)\Hx)  =  E [( Ac  (a,  fi)\ Hx  f  ]  - E2  [ Ac  (a,  fi)\ Hx \ , 
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=  Z Z  («JV  -  x> Pn  -  y)  S2  (a  -  X, P- -  y) 

x  y 

+  YjYjIM{^y)s2{oc-x,p- y) 

x  y 

Var( Ac  (aN,  PN  )|  H, )  =  KNg 3  («,v  -  x,  ^  -  y) 

jc  y  . 

(43) 

+1X4  (*,  y)  («iv  -  ^  Ar  -  y) 

x  y 

The  variance  of  the  Gaussian-based  LRT  given  H0  after  going  through  similar  steps  as  in 
the  Hx  case  is: 

Var(XG(a,p)\H^  =  YJYJ\jE[d(xi’y\)\81(a-x,p-y)\, 

x  y 

Var  (Ac  (aN,pN)\H0)  =  YJYJIM{^y)s1{aN-^PN-y)-  (44) 

x  y 

As  with  the  variance  of  the  Poisson-based  LRT  given  the  two  hypotheses,  the  variance  of 
the  Gaussian-based  LRT  given  both  hypotheses  are  approximately  equal  for  small  NEOs. 
The  SNR  of  the  Gaussian-based  LRT  is  bounded  by 

<82(aN-x,PN-y)  ZZLs!K"'i.A-j)  .(45) 

- - — -  - ■  <  SNRa  < 

{aN  ~y)  +  ^^IM{x,y)g2  (ciN  -x,/}N  -  y)  “  (x,y)g2  (aN -x,/3N  -  y) 

Inspection  of  Eq.  (40)  and  Eq.  (45)  indicates  that  the  two  algorithms  should  have  close  to 
the  same  level  of  performance  if  the  background  master  image  is  flat  and  small,  meaning 
that  no  background  stars  are  in  the  image.  This  allows  the  master  image  term  to  be  pulled 
out  of  the  summation  and  distributed  through  the  denominator,  which  results  in  the  same 
expression  for  SNR  whether  Gaussian  or  Poisson  distributions  are  used  for  photon 
statistics.  Both  expressions  of  the  bounds  on  the  SNRs  are  completely  characterized  by 
three  parameters:  the  number  of  photons  received  from  the  NEO,  the  background  master 
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image,  and  the  PSF.  The  next  section  will  introduce  models  used  to  simulate  these 
parameters  for  numerical  evaluation  of  the  SNR  expressions. 

3.6  Monte  Carlo  Simulations 

Simulations  developed  in  Matlab®  simulation  software  from  Mathworks®  will  be 
used  to  verify  the  feasibility  of  the  proposed  Poisson-based  LRT  compared  to  the 
Gaussian-based  LRT.  The  Matlab®  code  is  located  in  the  Appendix  and  is  well 
commented  to  provide  a  functional  description.  The  simulation  uses  the  most  detailed 
and  accurate  specifications  available  to  build  a  model  of  the  optical  system,  while  making 
appropriate  assumptions  where  necessary. 

The  code  first  creates  the  total  impulse  response,  to  include  a  model  of  PS  1  and 
the  average  effects  of  the  atmosphere.  The  code  then  generates  two  sets  of  data,  one  with 
a  simulated  NEO  present  and  one  without  a  simulated  NEO  present  using  Poisson 
statistics,  and  finally  it  calculates  the  LRT  based  on  Gaussian  assumptions  as  well  as  with 
Poisson.  The  impulse  response  is  used  to  create  simulated  images  by  convolving  it  with  a 
two  dimensional  Dirac  delta  function  scaled  by  the  expected  number  of  photons  for  the 
object  being  represented.  A  master  image  is  created  in  this  way,  but  simulated  data  is 
created  by  feeding  the  average  image  into  a  Poisson  noise  generator.  One  thousand  trials 
were  run  to  create  ROC  curves  for  each  of  the  two  algorithms  based  on  images  of  a  140 
meter  diameter  NEO  separated  from  a  star  by  28.21  arc-seconds.  The  angular  separation 
of  28.21  arc-seconds  equates  to  a  70  pixel  on  the  diagonal  separation  in  the  modeled 
detector  plain,  given  that  PS1  has  an  angular  resolution  of  0.285  arc-seconds  [11]. 
Different  ROC  curves  were  generated  for  star  magnitudes  ranging  from  0,  the  relative 
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brightness  of  Vega  to  25,  approximately  2.5  times  dimmer  than  a  140  meter  diameter 
NEO  at  1.3  AU  from  Earth.  A  simple  comparison  between  curves  can  provide  a  metric 
of  the  effectiveness  of  that  algorithm.  The  algorithm  with  a  ROC  curve  with  more  area 
under  the  curve  is  generally  a  better  perfonning  algorithm. 

Simulations  also  include  varying  other  environmental  situations  besides  the  size 
of  the  NEO.  The  magnitude  of  background  stars  in  the  image  is  varied,  as  well  as  the 
angular  separation  between  a  NEO  and  a  background  star  in  the  image.  Additionally,  ten 
thousand  trials  were  run  to  calculate  the  probability  of  false  alann  in  these  environments. 
The  increase  in  the  number  of  trials  for  false  alann  rates  provides  more  precision  in 
calculating  the  false  alarm  rate.  Limits  in  computing  power  prevent  Monte  Carlo 
simulations  as  a  method  to  calculate  a  threshold  that  would  yield  extremely  small  false 
alann  rates,  but  simulations  were  used  to  calculate  as  precise  a  rate  as  possible  with  the 
ten  thousand  trials.  False  alarm  rates  smaller  than  those  calculated  with  Monte  Carlo 
simulations  can  be  estimated  by  assuming  the  value  of  the  LRT  is  Gaussian  by  using  the 
Central  Limit  Theorem  [20],  Any  false  alann  rate  may  then  be  estimated  by  simply 
estimating  the  mean  and  variance  of  the  LRT.  Choosing  a  specific  false  alann  rate,  these 
three  plots  will  be  generated:  the  probability  of  detection  versus  angular  separation 
between  the  NEO  and  a  background  star;  the  probability  of  detection  versus  the 
magnitude  difference  between  the  NEO  and  a  background  star,  while  the  magnitude  of 
the  NEO  is  fixed;  and  the  probability  of  detection  versus  the  magnitude  of  the  NEO  while 
the  magnitude  of  the  background  star  is  fixed. 
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3.7  System  Model 

If  the  telescope  is  assumed  to  be  a  linear,  space-invariant  system,  all  that  is 
needed  to  model  the  telescope  is  its  impulse  response.  The  total  impulse  response  is  what 
is  imaged  when  the  system  is  interrogated  with  a  point  source.  The  system  includes  the 
atmosphere  and  the  telescope.  The  total  impulse  response  is  the  inverse  Fourier 
transform  of  the  product  of  the  average  atmospheric  transfer  function  with  the  optical 
transfer  function  of  the  telescope,  given  by 

*„(«.v)  =  F-‘{G„(/„/J,)G^(/„/,)} .  (46) 

The  first  thing  to  consider  when  building  the  system  model  is  proper  spatial 
sampling  for  the  detector  plane,  the  aperture  plane,  and  the  atmosphere.  If  not  sampled 
properly,  aliasing  will  cause  inaccurate  results  [16].  The  simplest  spatial  sampling  to 
detennine  is  in  the  detector  plane  due  to  the  fact  that  the  spatial  sampling  of  the  PS  1  in 
the  detector  plane  is  the  size  of  the  actual  pixels  of  the  CCD  array  of  PS1,  Ad  =10 jum 
[11]. 

The  spatial  sampling  in  the  aperture  is  slightly  more  difficult  to  arrive  at.  Nyquist 
sampling  theorem  requires  a  minimum  of  two  samples  per  period,  implying  that  the  phase 
change  in  the  aperture  plane  must  be  less  than  n .  The  field  at  the  aperture,  fa ,  due  to  a 
point  source  a  distance  z  away  from  the  aperture  is  defined  by 

j2m>[t-Ri(x,y)/c] 

^TT-TT - (47) 

jA,Rf(x,y) 
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where  x  and  y  are  points  measured  from  the  center  of  the  aperture.  R{  is  the  range 
from  the  point  source  to  the  aperture  as  a  function  of  the  location  in  the  aperture,  and  is 
described  by  the  distance  equation 

Ri  (*»  y)  =  t/K)2+(A)2+z2 .  (48) 


The  spatial  sampling  of  the  aperture,  Aa  is 


D 

~Ma 


D  is  the  diameter  of  the  aperture,  M a  is  the  number  of  samples  in  the  simulated 
aperture  plane.  The  binomial  approximation  of  Eq.  (48)  is 


Rl(x,y)~z 


(vA„)2  +(>’A,)2  ((»AT+(>A,)7 

2z2  Sz4 


(49) 


If  all  of  the  optics  of  PS1  are  treated  as  a  thin  lens,  the  phase  effects  of  the  field  through 
the  optics  of  the  telescope  can  be  modeled  at  the  aperture  with  a  lens  phase  screen 

t,enAX’y)  =  e  '  Af  (50) 


The  phase  of  the  lens  phase  screen  cancels  the  quadratic  term  in  Eq.  (49),  leaving  the 
third  tenn  as  the  largest  contributor  to  phase  change  at  the  aperture.  The  phase  change  of 
the  aperture  is  described  by 


n 


A  cp  =  ■ 


D 

V  2  J 


v  2 


4Af 


(51) 
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which  must  be  less  than  n  to  satisfy  Nyquist  sampling.  The  only  element  that  is  variable 
is  Afl  =  ,  specifically  Ma  must  be  large  enough  to  satisfy  the  phase  requirement. 


4T/3[l6  [  2M  a  )\ 


D4  J_ 
4  Af  16 


<1 


Using  500  mn  for  A  because  the  largest  phase  difference  will  occur  with  the  smallest 
wavelength,  8  m  for  the  focal  length  /  ,  and  1.8  m  for  the  diameter  D ,  Ma  >5123  was 


found  to  ensure  adequate  sampling  across  the  entire  aperture  [16]. 

The  spatial  sampling  of  the  atmosphere  requires  a  brief  description  of  the 
atmospheric  model.  The  largest  contributor  to  phase  error  due  to  a  turbulent  atmosphere 
is  tilt.  Tilt  can  be  modeled  as  a  zero  mean  Gaussian  random  variable,  which  is  correlated 
over  a  finite  amount  of  time.  Tilt  across  a  small  sample  of  time  is  highly  correlated,  but 
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samples  taken  over  longer  periods  of  time  become  more  and  more  uncorrelated,  until 
there  is  zero  correlation  between  samples.  The  exposure  time  of  one  image  taken  by  PS1 
is  typically  30  seconds.  Over  a  30  second  period,  the  atmosphere  at  the  start  of  the 
sample  time  is  uncorrelated  with  the  atmosphere  at  the  end  of  the  sample  time,  and  the 
effect  of  tilt  averages  to  zero.  This  allows  the  effect  of  tilt  to  be  ignored  in  this  model. 
Other  aberrations  caused  by  a  turbulent  atmosphere  can  be  modeled  by  an  average 
transfer  function.  The  average  transfer  function  for  a  long  exposure  is  given  by  [25] 


-3.44 

Gam{fxJy)  =  e 


^/A/X|/f+/y2 


(52) 


The  frequency  sampling  of  the  atmosphere  is 

A,  = — - — 

2  MA, 

where  M d  is  the  number  of  pixels  in  one  dimension  of  the  simulated  square  detector 
plane  and  was  chosen  to  be  100  pixels.  The  size  of  the  image  of  a  140  meter  wide  object 
at  1.3  AUs  with  a  focal  length  of  8  meters  on  the  detector  can  be  found  using  simple 
trigonometry  and  is  approximately  5.6  pm,  or  just  over  half  of  a  pixel;  therefore,  the 
complete  image  of  the  NEO  will  fit  well  within  100  pixels.  The  spatial  sampling  of  the 
atmosphere, 

a„=4/a/  =  ttt-  (53> 


ensures  that  Eq.  (52)  has  the  same  number  of  samples  in  the  frequency  spectrum  as  the 
transfer  function  of  the  telescope  [16]. 
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Now  that  the  proper  sampling,  Aatm,  Aa,  and  Ad  have  been  detennined  for  the 
model,  and  the  average  transfer  function  for  the  atmosphere  has  been  developed,  the  final 
step  in  obtaining  the  total  impulse  response  from  Eq.  (46)  is  to  develop  the  transfer 
function  of  the  telescope.  The  transfer  function  of  the  telescope  is  the  Fourier  transfonn 
of  the  PSF,  or  the  impulse  response  of  the  telescope.  This  is  detennined  by  propagating 
the  field  at  the  aperture  due  to  a  point  source,  through  the  optics  of  PS1.  First  a  circular 
array  that  represents  the  aperture  must  be  created.  The  field  at  the  aperture  can  be 
modeled  as  a  uniform  plane  wave.  This  is  an  appropriate  approximation  to  a  spherical 
wave  (point  source)  propagated  from  a  long  distance  [26],  If  the  field  at  the  aperture  is  a 
uniform  plane  wave,  then  the  aperture  may  be  modeled  as  a  binary  screen,  1  inside  the 
aperture  and  0  outside.  The  lens  phase  screen,  Eq.  (50),  is  then  multiplied  to  the  field  at 
the  aperture  to  create  the  field  f,  (xm,yn),  which  is  then  propagated  to  the  detector  plane, 
located  in  the  focal  plane  of  the  telescope  system.  The  Rayleigh-Sommerfeld  diffraction 
integral,  propagates  the  field  f,(xm,yn )  a  distance  /  ,  and  the  resulting  field  at  the 

detector,  fd(u,v,t )  is 


M, 


fd(u’v’  o=zz 


m= 1  n- 1 


jAR22(xm,yn,u,v) 


(54) 


in  which  R2  is  the  distance  between  every  pixel  in  the  aperture  plane  and  every  pixel  in 
the  detector  plane: 

Ri  (xm ,yn > v)  =  \l(xmAa  - uAd )2  +  ( ynAa  - vA df  +  f 2  (55) 
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The  squared  magnitude  of  Eq.  (54)  is  the  PSF  of  the  telescope,  so  [26] 

G,v(fxJr)  =  F^fd(uAd,vAd,tfY  (56) 

Now  Eq.  (46)  is  used  to  find  the  total  impulse  response  of  the  system  which  completely 
characterizes  the  system  for  an  arbitrary  input.  The  output  image  due  to  the  incidence  of 
any  input  can  be  modeled  as  the  convolution  of  the  input  with  the  impulse  response  [26]. 

3.8  Measured  Data 

Hardware  with  the  capabilities  of  PS1  is  not  available  to  test  the  proposed 
hypothesis  on  actual  image  data;  instead  a  relative  test  environment  is  developed  to  test 
the  difference  in  the  perfonnance  of  the  two  LRT  algorithms.  Having  a  smaller  telescope 
that  is  not  capable  to  detecting  NEOs,  a  brighter  object  is  used  and  the  integration  time  of 
the  telescope  will  be  adjusted  to  make  the  object  difficult  to  detect.  Polaris  is  chosen 
because  its  relative  fixed  location  in  the  sky  provides  time  to  collect  a  large  number  of 
images  without  having  to  readjust  the  hardware  or  to  register  the  images.  Also,  using  a 
bright  object  that  is  detectable  with  the  naked  eye  provides  a  known  truth  for  hypothesis 
testing.  Reducing  the  integration  time  of  the  CCD  device  on  the  telescope  reduces  the 
detectability  of  even  a  bright  object.  The  integration  time  will  be  set  to  such  a  level  that 
detection  is  difficult  and  a  statistically  significant  number  of  images  will  be  collected, 
ideally  a  minimum  of  one  thousand  images.  These  images  will  be  averaged  together  to 
verify  that  Polaris  is  detected  in  the  averaged  data.  Each  image  will  be  fed  into  each 
algorithm  and  ROC  curves  will  be  generated  and  compared  as  was  done  in  simulations. 
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Background  noise  images  from  a  location  in  the  image  that  does  not  have  an  object  will 
be  used  as  the  image  data  for  H0 . 

The  CCD  in  the  telescope  has  a  gain  factor  that  converts  the  photon  count  to  a 
digital  count  that  is  not  Poisson.  To  properly  implement  the  proposed  algorithm,  the 
digital  counts  must  be  converted  back  into  photon  counts.  For  this  conversion,  the  gain 
factor  of  the  CCD  is  required.  Consider  the  following. 

d  =(j)I 

where 

•  d  is  the  digital  count  data, 

•  (f)  is  the  CCD  gain  factor,  and 

•  1  is  the  Poisson  distributed  photon  count. 

The  expected  value  of  the  digital  count  is 

E[d]  =  <j>E[l],  (57) 

and  the  variance  of  the  digital  count  is 

Var(d)  =  </>2(E[l2~\-E2[lf) 

Var(d)  =  </>2Var(I)  =  t/>2E[l].  (58) 

The  ratio  of  (58)  to  (57)  results  in  an  expression  for  the  gain  factor, 

Var(d) 
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Both  the  variance  and  the  expected  value  of  the  digital  count  data  can  be  adequately 
estimated  with  a  reasonable  number  of  sample  images.  Converting  the  digital  count  back 
into  photon  count  will,  however  induce  some  quantization  noise  [24], 

3.9  Summary 

This  chapter  described  the  theory  supporting  the  basis  for  the  Poisson-based  LRT, 
as  well  as  the  methodology  for  comparing  the  perfonnance  of  the  Poisson-based  LRT  to 
the  perfonnance  of  the  Gaussian-based  LRT.  The  methodology  addresses  the 
comparison  of  the  two  algorithms  in  three  ways.  First,  analysis  is  performed  on  the 
expected  SNR  of  each  LRT.  The  LRT  which  has  a  larger  SNR  should  be  indicative  of 
the  better  perfonning  algorithm.  Secondly,  Monte  Carlo  simulations  are  performed  based 
on  complex  models  of  hardware  and  environmental  characteristics,  in  which  the  two 
algorithms  are  compared  in  several  different  scenarios.  Lastly,  measured  data  is  taken  in 
a  relative  test  environment  and  the  performance  of  each  algorithm  with  the  measured  data 
is  compared  using  ROC  curves.  The  results  of  these  methods  are  presented  in  the  next 
chapter. 
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4  Analysis  and  Results 


4.1  Chapter  Overview 

This  chapter  presents  the  results  of  the  three  approaches  described  in  Chapter  3  to 
comparing  the  Gaussian-based  LRT  to  the  Poisson-based  LRT.  The  first  of  these  results 
is  a  numerical  analysis  of  the  SNR  of  each  algorithm  based  on  either  actual  parameters  or 
modeled  parameters.  The  second  set  of  results  is  simulated  data  based  on  models  of  the 
hardware  and  environmental  characteristics.  The  third  set  of  results  is  measured  data 
taken  in  a  relative  test  environment. 

4.2  Analytic  Results 

The  SNR  analysis  requires  a  model  for  the  PSF,  a  model  for  the  expected  number 
of  photons  received  from  a  NEO  and  a  model  for  the  master  image.  The  PSF  was 
modeled  with  a  seeing  parameter  of  8  cm  as  described  in  the  previous  chapter,  and  details 
of  that  model  are  discussed  in  the  next  section  addressing  simulated  results.  The 
expected  number  of  photons  received  by  the  detector  of  PS  1  from  a  NEO  of  a  desired 
diameter  at  a  distance  of  1.3  AU  is  calculated  using  Eq.  (7)  from  Chapter  3.2.  The  master 
image  is  generated  by  convolving  point  sources  of  a  desired  magnitude  with  the  modeled 
PSF  and  adding  the  average  background  photons. 

The  upper  and  lower  bound  on  the  SNR  from  Eq.  (40)  and  Eq.  (45)  were 
calculated  using  modeled  inputs  for  various  master  images  and  varied  NEO  sizes.  The 
following  graphs  are  a  demonstrative  subset  of  the  results  for  the  varied  parameters. 
They  display  the  SNR  versus  the  apparent  magnitude  of  a  background  star  for  different 
angular  separations  between  the  modeled  NEO  and  star.  In  Figure  7  the  modeled 
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Figure  7.  SNR  of  140  m  NEO  vs.  Star  Mag; 


Figure  8.  SNR  of  140  m  NEO  vs.  Star  Mag; 


28.2  arc-sec  Angular  Separation 


14.1  arc-sec  Angular  Separation 


Figure  9.  SNR  of  400  m  NEO  vs.  Star  Mag; 


Figure  10.  SNR  of  400  m  NEO  vs.  Star  Mag;  0 


28.2  arc-sec  Angular  Separation 


arc-sec  Angular  Separation 


NEO  is  140  meters  in  diameter,  and  the  angular  separation  between  the  star  and  the  NEO 

is  approximately  28.21  arc-seconds.  The  star’s  apparent  magnitude  was  varied  from  0  to 

25.  Figure  8  shows  a  similar  set-up  where  the  only  change  is  the  angular  separation  to 

around  14.11  arc-seconds.  Figure  9  increases  the  size  of  the  NEO  to  400  meters  in 

diameter,  at  an  angular  separation  of  28.21  arc-seconds.  Since  the  SNR  for  such  a  large 

NEO  is  so  large  when  far  enough  away  from  a  nearby  star,  Figure  10  shows  the  predicted 

SNR  when  the  400  meter  diameter  NEO  is  in  the  same  pixel  as  a  background  star,  or  an 
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angular  separation  of  zero  arc-seconds.  These  results  predict  that  in  a  flat,  dim 
background  the  two  LRTs  should  perfonn  equally  as  well,  but  the  proposed  Poisson- 
based  LRT  should  outperfonn  the  Gaussian-based  LRT  when  there  is  a  bright  object 
nearby  and  the  background  is  not  flat.  Chapter  3.5  established  that  using  the  GLRT  could 
not  be  analyzed  using  SNR  calculations.  The  next  section  will  verify  the  results 
predicted  with  the  SNR  analysis  as  well  as  investigate  the  performance  of  the  GLRT 
through  simulations. 

4.3  Simulation  Results 

The  simulated  PSF  for  each  of  the  four  wavelengths  used  to  represent  Pan-Starrs’ 
wideband  filter  is  shown  in  Figure  11.  The  wavelength  dependent  PSFs  were  simulated 
using  a  r0  of  8  cm.  These  PSFs  were  used  to  simulate  images  for  testing  the  performance 
of  each  algorithm.  Simulated  images  were  created  by  convolving  the  four  PSFs  with  the 
object,  adding  the  four  resulting  images  together,  and  adding  Poisson  random  noise.  Both 
the  Gaussian  and  the  Poisson  LRT  require  the  use  of  the  PSF,  for  this  parameter  the 


(a)  (b) 


Figure  11.  Impulse  Response  for  (a)  500  nm,  (b)  600 
nm,  (c)  700nm,  and  (d)  800  nm  wavelengths 
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average  of  the  four  wavelength-dependent  PSFs  is  used.  In  all  Monte  Carlo  simulations, 
the  NEO  was  modeled  at  a  distance  of  1 .3  AUs 

Figure  12  shows  the  ROC  curves  generated  with  a  Magnitude  20  star,  28.21  arc- 
seconds  away  from  a  140  meter  NEO.  These  simulated  ROC  curves  support  the 
analytically  predicted  perfonnance  of  the  two  LRTs  with  a  dim  star  in  the  background, 
where  they  both  perform  equally  as  well.  Also  in  Figure  12  is  verification  that  for  dim 
background  stars,  the  perfonnance  of  the  Poisson-based  GLRT  where  the  unknown 
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Figure  12.  ROC  Curve;  Magnitude  20  Star 


Figure  13.  ROC  Curve;  Magnitude  10  Star 


Figure  14.  ROC  Curve;  Magnitude  0  star 
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parameter  KN  is  estimated  is  nearly  equivalent  to  the  other  two  LRTs.  If  the  apparent 
magnitude  of  the  star  in  the  background  is  decreased,  making  it  brighter,  the  performance 
of  the  Gaussian-based  LRT  should  degrade,  according  to  analytical  results.  Also,  the 
performance  of  the  Poisson-based  GLRT  should  degrade  due  to  the  increased  variance 
from  the  maximum  likelihood  estimate  as  discussed  in  Chapter  3.5.  Figure  13  and  Figure 
14  demonstrate  this  phenomenon  with  the  simulated  ROC  curves  for  the  same  NEO  at  the 
same  angular  separation,  but  with  a  Magnitude  10  and  Magnitude  0  star,  respectively. 

Figure  15  is  a  plot  of  the  detection  probability  versus  the  apparent  magnitude  of 
the  NEO.  The  NEO  magnitude  was  varied  from  21.79  (400  m  diameter)  to  24.8  (100  m 
diameter),  the  false  alarm  rate  was  set  to  lxl O'4,  and  the  detection  probability  was 
calculated  using  1000  simulated  realizations  of  images  with  Poisson  noise.  There  was  no 
star  simulated  in  the  scene,  only  background  radiation.  Figure  15  shows  there  is  nearly  a 
doubling  in  the  detection  probability  for  a  magnitude  24  NEO,  from  0.291  to  0.589. 
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Figure  15.  Detection  Probability  vs.  NEO  Magnitude 
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Difference  in  Magnitude  between  NEO  and  Star 

Figure  16.  Detection  Probability  vs.  NEO-Star  magnitude  difference 

Figure  16  is  a  plot  of  the  detection  probability  versus  the  difference  in  magnitude 
between  the  simulated  NEO  and  the  simulated  star  at  an  angular  distance  of  28.21.  The 
magnitude  of  the  NEO  was  held  constant  at  24,  while  the  magnitude  of  the  star  was 
varied  from  0  to  25.  Note  that  a  difference  of  5  in  apparent  magnitude  equates  to  a 
difference  in  relative  brightness  of  a  factor  of  100.  Figure  16  shows  that  the  detection 
probability  is  marginally  improved  by  the  proposed  algorithm  when  the  magnitude  of  the 
star  is  less  than  100  times  as  bright  as  the  NEO.  When  the  magnitude  of  the  star  is 
between  100  times  and  roughly  25,000  times  as  bright  as  the  NEO,  the  detection 
probability  is  approximately  doubled  by  the  proposed  algorithm.  When  the  magnitude  of 
the  star  is  greater  than  roughly  25,000  times  as  bright  as  the  NEO,  the  improvement  in 
detection  probability  by  the  proposed  algorithm  is  over  a  factor  of  4.  A  star  that  is  25,000 
times  as  bright  as  a  NEO  may  sound  like  an  unlikely  encounter,  but  this  is  only  a 
difference  in  magnitude  of  1 1,  and  in  this  example  the  star’s  magnitude  is  14.  There  are 
15.5  million  stars  as  bright  as  or  brighter  than  a  magnitude  14. 


Poisson  Model 
Gaussian  Model 
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Figure  17.  Detection  Probability  vs.  Angular  Separation 

The  final  simulated  environment  varied  the  angular  separation  of  a  200  meter 
diameter  NEO  and  a  magnitude  15  star  from  0  to  28.21  arc-seconds.  Again,  the  detection 
threshold  in  this  simulation  was  set  to  provide  a  false  alarm  rate  of  lxlCT4 .  Figure  17  is 
a  plot  of  detection  probability  versus  angular  separation  and  shows  that  a  NEO  of  this 
size  can  be  detected  with  a  probability  greater  than  0.9  approximately  3  arc-seconds 
closer  to  a  magnitude  15  star  with  the  proposed  Poisson-based  algorithm  than  with  the 
current  Gaussian-based  algorithm. 

4.4  Measured  Data  Results 

Two  sets  of  measured  data  were  collected  on  two  different  nights  with  two 
different  CCDs.  The  first  set  of  data  was  taken  in  early  August  in  Dayton,  OH.  The 
telescope  was  brought  into  focus,  and  the  viewfinder  was  calibrated  using  the  moon  as  a 
reference  due  to  the  ease  with  which  the  moon  is  located  with  an  un-calibrated 
viewfinder.  Polaris  was  then  located  and  its  identification  was  verified  by  its  relative 
location  to  the  group  of  stars  in  the  constellation  Ursa  Major,  known  as  the  Big  Dipper,  as 


57 


Figure  18.  Sample  image  of  Polaris  with  short  Figure  19.  Average  Image  of  Polaris 

exposure 

well  as  by  its  stationary  location  relative  to  the  rotation  of  the  Earth.  Once  Polaris  was 
lined  up,  347  images  were  taken  at  a  1  second  exposure  time.  These  long  exposures  were 
used  to  estimate  the  average  PSF  of  the  telescope  and  CCD.  The  integration  time  of  the 
CCD  was  then  reduced  to  10 jus,  the  point  at  which  Polaris  was  difficult  to  detect,  and 
253  images  were  taken  before  hardware  malfunctions  prevented  further  data  collection 
for  this  test.  The  gain  factor  was  calculated  by  subtracting  the  average  background  noise 
from  each  image  and  creating  an  average  “super-pixel”  around  the  location  of  Polaris  in 
the  image.  The  mean  and  variance  of  this  super-pixel  was  then  used  to  calculate  the 
estimated  gain  factor  of  0.7103,  indicating  that  this  particular  CCD  is  not  operating  in 
avalanche  mode,  and  it  takes  more  than  one  photon  to  produce  an  electron.  This  gain 
factor  was  removed  from  the  data.  Figure  18  shows  a  sample  image  taken  with  the 
integration  time  of  the  CCD  set  to  10  jus  .  Figure  19  shows  the  average  of  the  253  images 
taken  at  that  exposure  time,  and  indicates  the  presence  of  Polaris. 
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Figure  20.  ROC  Curve  for  Polaris  Short  Exposure  Images 
Figure  20  shows  the  ROC  curve  developed  by  using  the  10  /us  exposure  images  in 
the  two  algorithms.  As  with  the  simulated  data,  the  test  data  depicts  an  increase  in 
perfonnance  by  using  the  proposed  algorithm.  The  detection  probability  with  the 
Poisson-based  LRT  at  the  smallest  measureable  false  alarm  rate  is  0.5785,  while  the 
detection  probability  with  the  Gaussian-based  LRT  for  the  same  false  alarm  rate  is 
0.0807.  Based  on  this  limited  sample  set,  the  increase  in  detection  probability  for  very 
small  false  alarm  rates  is  by  a  factor  of  seven. 

The  second  set  of  data  was  taken  in  early  December  in  Dayton,  OH.  The 
telescope  was  again  brought  into  focus,  and  the  viewfinder  was  calibrated  using  the  moon 
as  a  reference  as  with  the  first  data  collection.  Immediately  it  was  apparent  that  the 
seeing  parameter  on  this  night  was  much  better  than  the  night  of  the  previous  data 
collection  as  many  more  stars  were  visible  that  were  not  visible  with  the  same  optics 
before.  A  star  in  the  general  location  of  where  Polaris  should  be  was  focused  on  that 
appeared  to  be  stationary  with  respect  to  the  rotation  of  the  Earth,  so  it  was  assumed  that 
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this  was  Polaris.  A  series  of  46  images  at  1  second  exposure  time  and  1179  images  at  10 
/us  exposure  times  were  taken.  Once  the  data  was  processed,  there  was  doubt  that  this 
star  was  Polaris  because  there  was  some  steady  movement  of  the  star.  This  required  that 
the  images  be  registered  before  they  could  be.  Another  issue  discovered  after  post¬ 
processing  the  data,  was  that  the  long  exposure  images  saturated  the  CCD,  this  means 
that  these  46  images  were  not  the  shape  of  the  PSF.  Instead  the  cross-section  of  the 
intensity  of  these  images  resembled  what  is  commonly  referred  to  as  a  “Top-Hat” 
function  where  the  top  portion  of  the  image  intensity  data  was  cut  off.  These  images 
were  used  to  estimate  the  PSF  by  creating  a  simulated  PSF  that  had  the  same  shape  and 
width  as  the  base  of  the  average  of  these  46  images.  In  a  similar  manner  as  with  the  first 
set  of  data,  the  estimated  gain  factor  was  found  to  be  1.1932,  indicating  that  this  CCD  is 
operating  in  avalanche  mode  and  more  than  one  electron  is  produced  for  every  photon. 
This  gain  factor  was  removed  from  the  data.  Figure  21  shows  a  sample  image  taken  with 
the  integration  time  of  the  CCD  set  to  10  jus  .  Figure  22  shows  the  average  of  the  1179 
images  taken  at  that  exposure  time,  and  indicates  the  presence  of  the  star. 


Figure  21.  Sample  Image  of  Star  with  Short 
Exposure 


Figure  22.  Average  Image  of  Star 
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Figure  23.  ROC  Curve  of  Second  Data  Set 

Figure  23  shows  the  ROC  curve  developed  by  using  the  10/w  exposure  images 
with  the  Poisson  LRT,  the  Poisson  GLRT,  and  the  Gaussian  LRT.  The  improvement  is 
not  as  significant  as  with  the  first  set  of  data  and  there  are  several  possible  contributing 
factors  for  this.  The  first  set  of  data  had  a  limited  number  of  samples;  whereas  the  second 
set  of  data  had  a  much  more  statistically  significantly  number  of  samples,  also  r0  on  the 
second  night  of  collection  was  much  better,  making  the  average  PSF  smaller.  The  second 
set  of  data  with  more  samples  more  closely  resembles  the  simulated  performance  of  a 
NEO  without  the  presence  of  a  star  nearby.  Analysis  and  simulations  done  with  models 
of  Pan-Starrs  predicted  the  performance  of  the  LRTs  and  the  GLRT  should  be 
approximately  the  same  given  that  Polaris  was  the  only  thing  in  the  image.  The  detection 
probability  with  the  Poisson  LRT  at  the  smallest  measureable  false  alarm  rate  is  0.9101, 
while  the  detection  probability  with  the  Gaussian  LRT  for  the  same  false  alarm  rate  is 
0.81 17.  The  nearly  10%  increase  in  probability  of  detection  realized  in  this  data  set  could 
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indicate  the  resilience  of  the  proposed  Poisson-based  LRT  to  other  sources  of  noise  not 
accounted  for  in  simulations  and  SNR  analysis.  One  such  source  of  noise  is  scintillation 
in  which  the  PSF  fluctuates  over  time.  The  10  jus  integration  time  for  the  measured  data 
would  see  a  larger  fluctuation  in  the  shape  of  the  PSF  than  the  30  seconds  modeled  for 
Pan-Starrs  simulations,  the  PSF  of  which  would  more  closely  resemble  the  average  PSF. 
Further  research  is  required  to  definitively  characterize  the  realized  improvement  in  this 
data  set  versus  the  expected  performance  based  on  simulations  and  analysis. 

4.5  Summary 

This  chapter  presented  the  numerical  and  graphical  results  of  the  analysis, 
simulation,  and  measurement  of  the  performance  of  the  proposed  algorithm  in 
comparison  with  the  existing  matched  filter,  or  Gaussian-based  algorithm.  The  results  in 
all  cases  studied  demonstrate  an  improved  performance  by  the  proposed  algorithm  over 
the  matched  filter. 
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5  Conclusions  and  Recommendations 


5.1  Chapter  Overview 

This  chapter  presents  conclusions  and  recommendations  derived  from  perfonning 
this  research  effort.  The  chapter  begins  by  summarizing  conclusions  drawn  from  the 
results  presented  in  Chapter  4.  Next,  the  significance  and  impact  of  the  inferred 
conclusions  are  projected.  Finally,  recommendations  are  made  for  immediate  action  and 
for  future  work. 

5.2  Conclusions  of  Research 

The  Poisson-based  LRT  algorithm  proposed  by  this  research  produced  increases 
in  probability  of  detection  as  high  as  a  factor  of  seven  over  the  existing  algorithm  for 
measured  data.  In  all  simulated  conditions  explored,  the  proposed  algorithm  perfonned 
as  well  as  the  current  matched  filter  algorithm.  In  certain  simulated  conditions  such  as  a 
dim  NEO  near  a  bright  star,  the  detection  rate  for  small  false  alann  rates  was  more  than  a 
factor  of  four.  Little  has  changed  in  the  past  few  decades  with  regard  to  the  basic  image 
signal  processing  theory  used  in  the  detection  of  dim  astronomical  objects.  The  optical 
systems  have  increasingly  become  more  capable  and  functional,  but  how  the  data 
received  from  the  optical  systems  are  processed  remains  fundamentally  the  same,  with 
the  exception  of  the  Pan-Starrs  program  and  the  projected  LSST  program. 

5.3  Significance  of  Research 

The  significance  of  this  research  is  potentially  very  important  to  the  astronomical 

community.  The  research  demonstrated  improvement  in  binary  detection  applications 

through  theoretical  analysis,  simulated  hardware  used  in  NEO  detection,  and  with  real 
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images  taken  with  two  different  cameras  on  a  much  smaller  telescope  than  would  feasibly 
be  used  for  NEO  detection.  This  leads  to  the  conclusion  that  the  proposed  algorithm  has 
utility  across  all  applications  of  detecting  the  presence  of  any  dim  astronomical  object 
using  terrestrial-based  optical  systems.  Potential  significance  of  the  research  ranges  from 
an  increased  capability  of  tracking  debris  in  orbit,  to  detecting  a  previously  undetectable 
asteroid  or  comet  on  a  collision  course  with  Earth  early  enough  to  establish  an  effective 
plan  of  action  to  save  millions  of  lives. 

5.4  Recommendations  for  Action  and  Future  Research 

The  research  conducted  here  was  unfunded;  however,  based  on  preliminary 
results  prior  to  the  completion  of  this  research,  interest  from  one  particular  survey 
program  was  expressed  in  the  potential  impacts  of  this  research.  There  are  a  few 
challenges  remaining  before  the  proposed  algorithm  is  easily  incorporated  to  existing 
hardware.  A  thorough  study  of  calibration  issues,  hardware  specific  parameters,  and  any 
other  engineering  challenges  in  order  to  optimize  the  algorithm  is  required.  A  method  is 
needed  that  determines  the  threshold  of  detection,  which  varies  from  one  location  in  the 
sky  to  another  due  to  the  background  image,  and  is  also  hardware  specific.  Once  these 
challenges  are  overcome,  the  research  should  be  extended  to  measure  the  actual  gain 
realized  on  more  sophisticated  astronomical  imaging  equipment. 

5.5  Summary 

In  the  final  chapter,  conclusions  drawn  from  the  research  conducted  are  presented. 
Potential  significance  of  the  results  reported  is  forecasted  and  recommendations  for 
actions  based  on  this  research  as  well  as  future  research  to  be  conducted  are  presented. 
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Appendix 


Al.  Simulation  MatLab  Code 

%Capt  Curtis  Peterson  MatLab  M-Code 
%Pan-Starrs  model  and  Simulations 

%%  Adjustable  Parameters 

Target_radius  =  70; 
trials  =  1000; 
rho_t  =  .03; 
r_0  =  0.08; 

Star_Mag  =  10; 

NEO_position  =  [15  15]; 
star_position  =  [85  85]; 
masterjmages  =  100; 
run_psl_psf  =  0; 


estimate_K_N  =  1; 

%%  Fixed  Parameters 
M_a  =  5123; 

M_d  =  100; 

D  =  1.8; 

f  =  8; 

dx  =  D/M_a; 
dy  =  dx; 
dxx  =  10e-6; 
dyy  =  dxx; 

AU  =  1.495978707ell; 

Target_Range  =  1.3*AU; 
h  =  6.626e-34; 
deltajam  =  .1; 
delta_t  =30; 

K_B  =  7; 
rl  =  M_a/2; 
r2  =  0; 

total_imp_resp  =  zeros(M_d,M_d,4];  %  allocate  space  for  impulse  response 

psl_psf  =  total_imp_resp;  %  allocate  space  for  psf 

K_N  =  zeros(l,4];  %  allocate  space  for  expected  number  of 

%  photo-electrons  received  from  NEO 

vega  =  [5.53215ell  5.1066ell  4.72506ell  4.38756ell];  %photons  received  by 

%1.8  m  aperture  from 
%Vega  per  wavelength 

K_S  =  vega/(2.512AStar_Mag];  %photo-electrons  received  front  nearby  star  per 

%wavelength 

angluar_separation  =  norm(star_position-NEO_position]*0.285;  %starto  NEO 

%separation  in  arc-seconds 

for  wavelength  =  1:4; 

%%  Create  Impulse  Response 

%  select  the  wavelength  dependent  parameters 
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%  NEO  radius  in  meters 

%  number  of  trials  for  Monte  Carlo  Simulations 
%  Target  reflectivity 

%  atmospheric  seeing  parameter  in  meters 
%  Apparent  Magnitude  [AM]  in  AB  system  of  nearby  star 
%  Pixel  location  (x,y)  of  the  NEO 
%  Pixel  location  (x,y)  of  nearby  star 
%number  of  images  to  average  for  the  master  image 
%  1  ==  run  code  to  develop  psf  of  psl  (very  time 
%  consuming],  0  ==  skip  this  section  of  code  (if  you 
%  have  previously  run  it  and  saved  the  psf.mat] 

%  1  ==  run  code  to  estimate  the  K_N  parameter  in  the 
%  Poisson  LRT,  0  ==  set  K_N  equal  to  1. 

%  number  of  samples  in  aperture 
%  number  of  samples  (pixels]  in  the  simulated  detector 
%  diameter  of  the  mirror  (aperture]  in  meters 
%  focal  length  of  the  system  in  meters 
%  sample  size  in  the  aperture  in  meters 

%  sample  size  in  the  detector  in  meters 

%  1 AU  in  meters 

%  range  to  target  in  meters  (definition  of  NEO] 

%  Planck's  constant 

%  Bandwidth  of  receiver  in  units  of  micrometers 
%  exposure  time  of  PS1  in  seconds 
%  average  background  photo-electrons  per  pixel 


if  wavelength  ==  1 
lam  =  0.5e-6; 

S_irr  =  2e3; 
tau_opt  =  0.67; 
tau_atm  =  .98; 
elseif  wavelength  ==  2 
lam  =  0.6e-6; 

S_irr  =  1.8e3; 
tau_opt  =  0.85; 
tau_atm  =  .99; 
elseif  wavelength  ==  3 
lam  =  0.7e-6; 

S_irr  =  1.4e3; 
tau_opt  =  0.99; 
tau_atm  =  .9; 
else  lam  =  0.8e-6; 

S  jrr  =  le3; 
tau_opt  =  0.95; 
tau_atm  =  .9; 
end 

if  run_psl_psf==l 

%create  M_a  by  M_a  aperture  array  with  a  circular  binary  screen  equal  to  the 
%size  of  the  mirror 
mi  =  floor(rl)  +  1; 
aperture_array  =  zeros(M_a,M_a); 
for  ii  =  l:M_a 
for  jj  =  l:M_a 

dist  =  sqrt((ii-mi)A2+(jj-mi)A2); 
if(dist  <=  rl) 
if(dist  >=  r2) 

aperture_array(ii,jj)=l; 

end 

end 

end 

end 

x  =  dx*(-floor(rl):floor(rl)); 

xx_mat  =  ones(M_a,l)*x;  %create  matrices  for  vectorizing  the 

%lens  phase  and  field  propagation 
%math  in  the  next  two  steps 

yy_mat  =  x'*ones(l,M_a); 

%Treat  the  telescope  as  a  single  thin  lens  with  focal  length  f  and 
%apply  the  lens  phase  to  the  aperture 
lens_phase  =  -pi*(xx_mat.A2  +  yy_mat.A2)/(Plam); 
source_array  =  aperture_array.*exp(lj.*lens_phase]; 

%Propagate  the  aperture  field  to  the  receiver  array  using  the 
%Rayleigh-Sommerfeld  diffraction  integral 
receive  r_array  =  zeros  (M_d,M_d); 
for  xx  =  l:M_d 

xxc  =  (xx  -  ceil(M_d/2])*dxx; 
foryy  =  l:M_d 

yyc  =  (yy  -  ceil(M_d/2])*dyy; 

R  =  (fA2  +  (xx_mat-xxc).A2  +  [yy_mat-yyc).A2).A(0.5]; 


%wavelength  in  meters 
%solar  spectral  irradiance  (W/m/um) 
%optical  transmittance  of  psl 
%atmospheric  optical  transmittance 
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receiver_array(yy,xx)  =  sum(sum(dx*dy*source_array.*... 
exp(2*pi*lj.*R./lam)))./(lam*lj*f); 

end 

end 

%find  the  point  spread  function  of  PS1,  which  is  normalized 
psl_psf(:,:, wavelength)  =  abs(receiver_array).A2; 
psl_psf(:,:, wavelength)  =  psl_psf(:,:, wavelength)/... 

sum(sum(psl_psf(:,:, wavelength))); 

end 

%  Find  the  average  atmospheric  transfer  function  for  a  long  exposure 
%  and  multiply  it  by  the  transfer  function  of  PS1  in  order  to  get  the 
%  total  transfer  function.  Take  the  inverse  Fourier  transform  of  the 
%  total  transfer  function  in  order  to  find  the  total  impulse  response 
dx_otf  =  lam*f/ (M_d*dxx); 
avg_otf  =  zeros(M_d,M_d); 
mii  =  floor(M_d/2)  +  l; 
for  i  =  l:M_d 
for  j  =  l:M_d 

dist  =  sqrt((i-mii)A2+(j-mii)A2); 
if(dist<=2*(D/2)/dx_otf) 

avg_otf(i,j)=exp(-3.44*((dist/(r_0/(dx_otf)))A(5/3))); 

end 

end 

end 

tot_otf  =  fftshift(avg_otf).*(fft2(fftshift(psl_psf(:,:, wavelength)))); 
totalJmp_resp(:,:, wavelength)  =  abs(ifftshift(ifft2(tot_otf))); 

%%  Calculate  Expected  number  of  photons  from  NEO  for  each  wavelength 

v=3e8/lam;  %frequency 

%dA  is  the  smallest  of  the  IFOV  and  the  area  of  the  target 

1F0V  =  (dxx/PTarget_Range)A2; 

A_B  =  Target_radiusA2*pi; 
if  IFOV  <  A_B 
dA  =  IFOV; 
else 

dA  =  A_B; 
end 

K_N  (wavelength)  =  SJrr*delta_lam*dA*rho_t*tau_atm*tau_opt*DA2*... 
delta_t/ (4*T  arget_Range  A  2*h*v); 

end 

%%  SNR  Calculations 

%the  following  variables  and  expressions  occur  in  the  SNR  calculations 
psf_avg  =  mean(total_imp_resp,3);  %  average,  non-wavelength  dependent  psf 
R_h  =  max(max(xcorr2(psf_avg)));  %  zero  lag  auto  correlation  of  psf 
KN  =  sum(K_N);  %total  expected  photo-electrons  from  NEO 
!M_h2  =  sum(sum(I_M.*(circshift(psf_avg,NEO_position-50).A2))); 
h2_IM  =  sum(sum((circshift(psf_avg,NEO_position-50).A2)./l_M)); 

KN_h3  =  sum(sum(KN*psf_avg.A3)); 

h3_IM  =  sum(sum((KN*circshift(psf_avg,NEO_position-50).A3)./(I_M.A2))); 
Var_Gaus_Hl  =  KN_h3  +  IM_h2; 

Var_Gaus_HO  =  !M_h2; 

Var_Pois_Hl  =  h3_IM  +  h2_IM; 
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Var_Pois_HO  =  h2_IM; 

SNR_Gaus_upper(NEOsize,starsize]  =  (KN*RJi)/sqrt(Var_Gaus_HO); 
SNR_Gaus_lower(NEOsize,starsize)  =  (KN*Rji)/sqrt(Var_Gaus_Hl); 
SNR_Pois_upper(NEOsize,starsize)  =  (KN*h2JM)/sqrt(Var_Pois_H0); 
SNR_Poisjower(NEOsize,starsize)  =  (KN*h2JM)/sqrt(Var_Pois_Hl); 

%%  Find  Probabilities  of  detection  and  Probabilities  of  false  alarm  for  the 
%given  NEO  Magnitude,  Star  Magnitude,  Angular  Separation 

NEO_mag  =  log(sum(vega)/sum(K_N))/log(2.512);  %Apparent  Magnitude  of  NEO 
NEO  =  zeros(M_d,M_d,4); 

NE0(NE0_position(l),NE0_position(2),:)  =  K_N; 

%conv2  is  used  below  to  create  the  image  of  the  NEO  (similarly  to  create  the 
%image  of  the  star)  though  the  PS1  optics  as  opposed  to  using  fft2  due  to  the 
%periodicity  of  the  fft2  function  which  caused  edge  effects  and  inaccurate 
%results  when  either  the  star  or  the  NEO  was  placed  close  to  the  edge  of  the 
%simulated  image 

NEOJmage  =  conv2(NEO(:,:,l),totalJmp_resp(:,:,l),'same')  +  ... 
conv2(NEO(:,:,2),totaljmp_resp(:,:,2),'same')  +  ... 
conv2(NEO(:,:,3),total_imp_resp(:,:,3),'same')  +  ... 
conv2(NEO(:,:,4),total_imp_resp(:,:,4),'same'); 

[col  row]  =  find(NEO_image==max(max(NEO_image))); 

if  estimate_K_N  ==  0; 

eta_pois  =  0:le-4:3;  %Poisson  LRT  range  of  thresholds 
elseif  estimate_K_N  ==  1; 

eta_pois  =  [0:.001:10  10:. 1:100  100:1000  Ie3:10:le4  Ie4:le2:le5  ... 
Ie5:le3:le6  Ie6:le4:le7]; 

end 

eta_gaus  =  -100+200/length(eta_pois):... 

200/length(eta_pois):100;  %Gaussian  LRT  range  of  thresholds 
det_pois  =  zeros(length(eta_pois), trials); 
fa_pois  =  det_pois; 
fa_gaus  =  det_pois; 
det_gaus  =  det_pois; 

star  =  zeros(size(total  jmp_resp)); 
star(star_position(l),star_position(2),:)  =  K_S; 

starjmage  =  conv2(star(:,:,l),total_imp_resp(:,:,l),'same')  +  ... 
conv2(star(:,:,2),total_imp_resp(:,:,2),'same')  +  ... 
conv2(star(:,:,3),total_imp_resp(:,:,3),'same')  +  ... 
conv2(star(:,:,4),total_imp_resp(:,:,4),'same'); 

1_M  =  (starjmage  +  K_B); 


%the  variable  fx  below  is  a  fixed  value  for  all  trials  and  is  placed  outside  of  the  for  loop 
%of  trials  below  in  order  to  reduce  run  time 
fx=  circshift(psf_avg, [col-50  row-50]). /I_M; 

for  trial  =  l:trials; 
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d_0  =  poissrnd(I_M);  %  noisy  image  data  with  HO 

d_l  =  poissrnd(I_M  +  NEO_image);  %noisy  image  data  with  HI 
if  estimate_K_N  ==  1; 

K_N_hat_H0  =  sum(sum((d_0-I_M).*((d_0-I_M)>0))); 

K_N_hat_Hl  =  sum(sum((d_l-I_M).*((d_l-I_M)>0))); 
elseif  estimate_K_N  ==  0; 

K_N_hat_H0  =  1; 

K_N_hat_Hl  =  1; 
end 

%calculate  the  Poisson  LRT  values  for  both  hypotheses 

Poisson_LRT_Hl  =  sum(sum((d_l/  K_N_hat_Hl)  .*  log(l  +  K_N_hat_Hl*fx))) 
Poisson_LRT_H0  =  sum(sum((d_0  /  K_N_hat_H0).*  log(l  +  K_N_hat_H0*fx))) 
%calculate  the  Gaussian  LRT  values  for  both  hypotheses 
Gaussian_LRT_HO  =  conv2(d_0-I_M,psf_avg,'same'); 

Gaussian_LRT_HO  =  Gaussian_LRT_HO(col-l, row-1); 

Gaussian_LRT_Hl  =  conv2(d_l-I_M,psf_avg,'same'); 

Gaussian_LRT_Hl  =  Gaussian_LRT_Hl(col-l, row-1); 

%determine  false  alarms  and  detections  for  the  range  of  thresholds 
for  k  =  l:length(eta_pois) 

fa_gaus(k, trial)  =  Gaussian_LRT_HO  >=  eta_gaus(k); 
det_gaus(k, trial)  =  Gaussian_LRT_Hl  >=  eta_gaus(k); 
fa_pois(k, trial)  =  Poisson_LRT_H0  >=  eta_pois(k); 
det_pois(k, trial)  =  Poisson_LRT_Hl  >=  eta_pois(k); 
end 
end 

%Calculate  the  Probabilities  of  False  Alarm  and  Detection  for  the  given 

%number  of  trials 

Pd_Pois  =  sum(det_pois,2)/trials; 

Pd_Gaus  =  sum  (det_gaus, 2) /trials; 

Pfa_Pois  =  sum(fa_pois,2)/trials; 

Pfa_Gaus  =  sum(fa_gaus,2)/trials; 
exit 
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A2.  Measured  Data  MatLab  Code 


%Capt  Curtis  Peterson 

%Thesis  Measured  Data  MatLab  M-code 

%Compare  Poisson  Based  LRT  to  Gaussian  Based  LRT  with  images  of  Polaris 
%taken  at  a  lOus  exposure  time. 

estimate_K_N  =  1;  %choose  whether  to  estimate  the  value  of  K_N 
%(l==estimate) 

%%  Import  .jpg  image  data  (modify  for  other  image  formats) 

%Create  NxMxdummy  array  "A"  which  is  an  array  of  registered  and  %resolved  two  dimensional 
intensity  images  of  Polaris  from  a  long  %integration  time.  This  will  be  used  to  estimate  the  PSF  of  the 
%telescope. 

psfjpegs  =  dir(['<///ename>','*.jpg  ]); 

[N,  M,  dummy]  =size(imread(psfjpegs(l).  name)); 
psf_numfiles  =  length  (psfjpegs); 

A  =  zeros(N,M,psf_numfiles); 

for  k  =  l:psf_numfiles 
1  =  imread(psfjpegs(k).name); 

I  =  mean(l,3); 

A(:,:,k)  =  mat2gray(I); 

end 

%Create  NlxMlxdummyl  array  "images"  which  is  an  array  of  registered  %unresolved  two 
dimensional  intensity  images  of  Polaris  from  a  short  %integration  time.  This  will  be  used  as  a  noisy 
data  image  for  LRT  %detection 
imagesjpegs  =  dir(['</i7ename>','*.jpg']); 

[Nl,Ml,dummyl]=size(imread(images  jpegs(l).name)); 

images_numfiles  =  length(imagesjpegs); 
images  =  zeros(N,M,images_numfiles); 

fork=  l:images_numfiles 

II  =  imread(imagesjpegs(k).name); 
images(:,:,k)=mean(Il,3); 

end 

%%  Estimate  the  telescope  PSF 

%find  the  max  of  the  average  PSF  and  shift  it  to  the  zero  lag,  or  the 
%center  of  the  image,  then  create  a  40x40  image  of  the  PSF  less  the  average 
%background  noise,  and  normalize  it,  this  is  the  estimated  PSF.  A  40x40 
%array  is  used  because  that  is  the  size  of  the  region  of  the  image  that 
%will  be  searched  for  an  "object"  (Polaris)  from  the  noisy  data  images 
average_psf  =  mean(A,3);  %the  average  resolved  image  of  Polaris  has  the 
%general  shape  of  the  PSF 
m=max(max(average_psf)); 

[col,row]=find(average_psf==m); 
shiftR  =  N/2  -  col  +  1; 
shiftD  =  M/2  -  row  +  1; 

avg_noise  =  mean2(average_psf(:,l;row-100)); 
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psf_temp  =  circshift(average_psf,[shiftR,shiftD]); 
psf_test  =  psf_temp(N /2-20:N/2+ 19, M/2-20:  M/2+19); 
psf_test  =  (psf_test  -  avg_noise).*((psf_test  -  avg_noise]>0); 
psf_test  =  psf_test/max(max(psf_test)); 

%find  the  location  of  Polaris  in  the  noisy  data  set  by  averaging  them 
%together  and  finding  the  peak.  A  40x40  region  around  the  location  of 
%Polaris  will  be  used  as  the  search  region 
averagejmage  =  mean(images,3); 
ml=max(max(average_image)); 

[coll,rowl]=find(average_image==ml); 

shift_a  =  coll  -  20:coll  +  19; 
shift_b  =  rowl  -  20:rowl  +  19; 

%%  Calculate  the  CCD  Gain  factor 

%  the  gain  factor  is  estimated  as  the  average  gain  over  the  40  pixels 
%  of  interest  which  is  the  mean  over  all  the  images  of  the  average 
%  value  of  the  40x40  pixels  minus  the  mean  background  noise  for  the 
%  same  pixels  divided  by  the  variance  of  the  same 
imagesl=mean(mean(images(shift_a,shift_b,:],l),2]; 
for  k  =  l:images_numfiles 
noise(k)=mean2(images(:,l:rowl  - 100, k)); 
end 

images2(l,:)  =  imagesl(l,l,:); 
images2  =  images2  -  noise; 

Gain  =  mean(images2)./var(images2]; 
images  =  images.  *Gain; 

%%  Calculate  the  Poisson  and  Gaussian  LRTs  for  HI  and  HO 

Poisson_LRT_Hl  =  zeros(length(shift_a),length(shiftJ)),size(images,3)); 

Poisson_LRT_H0  =  Poisson_LRT_Hl; 

Gaussian_LRT_Hl  =  zeros(l,size(images,3)); 

Gaussian_LRT_H0  =  Gaussian_LRT_Hl; 

Poisson_LRT_Hl_Max  =  Gaussian_LRT_Hl; 

Poisson_LRT_H0_Max  =  Gaussian_LRT_Hl; 

for  im_num  =  l:size(images,3] 

d_l  =  images  (shift_a,shift_b,im_num);  %data  for  HI  for  image  number 

%im_num 

d_0  =  images (shift_a  +  150,shift_b  -  150,im_num);  %data  for  HO  for 

%image  number  im_num,  using 
%background  noise  from  some 
%offset  from  the  location  of  %Polaris 
I_M  =  mean2(d_0)*ones(size(d_0));  %the  "master  Image"  is  a  constant 
%equal  to  the  average  background 
%noise 

if  estimate_K_N  ==  1; 

K_N_hat_H0  =  sum(sum((d_0-I_M).’|!((d_0-l_M)>0))); 

K_N_hat_Hl  =  sum(sum((d_l-I_M).*((d_l-l_M)>0)]]; 
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elseif  estimate_K_N  ==  0; 

K_N_hat_H0  =  1; 

K_N_hat_Hl  =  1; 
end 

for  alpha  =  l:length(shift_a); 
for  beta  =  l:length(shift_b); 

%shift  the  PSF  through  every  pixel  in  the  image, 

%calculate  the  LRT  and  then  find  the  max  value  of  the 
%LRT  and  compare  it  to  a  threshold 
f  =  (circshift(psf_test,... 

[shift_a(alpha)-20  shift_b(beta)-20]]]./I_M; 

Poisson_LRT_Hl  (alpha,  beta,  im_num)  =  sum(sum((d_l/K_N_hat_Hl).*... 
log(l+K_N_hat_Hl*f))); 

Poisson_LRT_H0  (alpha,  beta,  im_num)  =  sum(sum((d_0/K_N_hat_H0).*... 
log(l+K_N_hat_H0*f])); 

end 

end 

%max  Pois  LRT  value  for  HI  data: 

Poisson_LRT_Hl_Max(im_num)  =max(max(Poisson_LRT_Hl(:,:,im_num))); 

%max  Pois  LRT  value  for  HO  data,  and  estimating  G: 
Poisson_LRT_HO_Max(im_num]  =  max(max(Poisson_LRT_HO(:,:,im_num)]); 

%max  Gaus  LRT  value  for  HI  data: 

Gaussian_LRT_Hl(im_num)=max(max(conv2((d_l-I_M),psf_test,'same')]); 


%max  Gaus  LRT  value  for  HO  data: 

Gaussian_LRT_H0(im_num)=max(max(conv2((d_0-I_M),psf_test,'same'))); 

end 

eta_gaus  =  0:. 1:6000;  %range  of  threshold  for  Gaussian  LRT 
if  estimate_K_N  ==  0; 

eta_pois  =  240:140/length(eta_gaus):380-140/length(eta_gaus); 

%Poisson  LRT  range  of  %thresholds 


elseif  estimate_K_N  ==  1; 

eta_pois  =  1.2e5:3e5/length(eta_gaus):4.51e5-3e5/length(eta_gaus); 

%range  of  threshold  for  %Poisson  LRT 
with  K_N  %estimated 


end 


%allocate  space  for  probabilities  of  detection  &  false  alarms  for  all  %tests 

pd  =  zeros(l,length(eta_gaus]]; 

pfa  =  pd; 

pd_g  =  pd; 

pfa_g  =  pd; 

%compare  peak  detection  of  the  tests  to  a  varying  threshold  to  build  a 
%ROC  curve 

fork=  1:  length  (eta_gaus) 

pd(k)=sum(Poisson_LRT_Hl_Max>eta_pois(k))/length(Poisson_LRT_Hl_Max); 
pfa(k)=sum(Poisson_LRT_HO_Max>eta_pois(k))/length(Poisson_LRT_HO_Max); 
pd_g(k)  =  sum(Gaussian_LRT_Hl>eta_gaus(k))/length(Gaussian_LRT_Hl); 
pfa_g(k)  =  sum(Gaussian_LRT_HO>eta_gaus(k))/length(Gaussian_LRT_HO); 
end 
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